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A Solution 
Equations for a Compressible Fluid with 
Variable Properties, Including Dissociation’ 


of the Laminar Boundary-Layer 


L. L. MOORET 
The RAND Corporation 


SUMMARY 


The laminar boundary-layer equations are solved for a flat 
plate in compressible flow where equilibrium dissociation of the 
air is assumed, as well as variable air properties. Evaluation 
of the skin-friction and heat-transfer characteristics is made for 
a Mach Number range | to 20. Both the insulated plate with- 
out radiation and the plate with heat transfer are considered. 

The nonlinear partial differential equations having been re- 
duced to five first-order ordinary equations in integral form are 
solved by the differential analyzer, which is ideal for handling 
nonlinear relations, as well as having the facility for absorbing 
graphical air property data. 

Values of the physical and thermodynamic properties of air 
expressed in the quantities Pr (Prandtl Number), pyu/ppu, 
(density-viscosity parameter), and j (enthalpy parameter) are 
obtained from experimental data for lower temperatures and with 
the aid of the kinetic theory at the higher temperatures. Equilib- 
rium dissociation produces an appreciable change in composition 
of air at the higher temperatures, thus introducing into enthalpy 
the heat of formation of dissociated products, as well as altering 
the viscosity, Prandtl Number, etc. 

Values of local skin-friction coefficient, cy, and heat-transfer 
coefficient, h., are computed for the plate with heat transfer from 
the solutions of the boundary-layer equations assuming dissoci- 
Insulated plate temperature, 7;, and skin-friction co- 
efficient for the insulated plate are also computed in a similar 


ation. 


Comparisons are made of these quantities with the 
The comparisons 


manner. 
values computed by others for no dissociation. 
show that, for heat transfer to the plate, cy and heat transmis- 
sion per unit area Q[= h.(7; T..)] are essentially unaffected 
by dissociation for plate temperatures in the range 400° to 2,500° 
R. The friction coefficient, c;, for the insulated plate is also un- 
affected by dissociation. However, there is an apparent de- 
crease in insulated plate temperature, 7;, with a corresponding 
increase in h, due to dissociation, both varying with ambient 


Received December 14, 1951. 

* This paper is based upon a research problem conducted by 
the author in September, 1948, for the Masters thesis in the field 
of Engineering at the University of California at Los Angeles. 

t Research Engineer, Missiles Division. 


pressure. For extremely high plate temperatures (i.e., 7. = 25), 
dissociation affects Q greatly and the skin friction to a lesser ex- 
All computations are based on a temperature outside of 
100°R. 
dissociation affecting the skin 

parameters is shown to be the presence of high temperature at 


tent. 
the boundary layer, 75 = In conclusion, the criterion for 


friction and heat-transmission 
the plate which causes dissociation at that location, although 
the skin friction for the insulated plate is an exception to this rule 


SYMBOLS 


velocity of sound,a = V ygRT, (ft. per sec.) 
local friction coefficient, ¢, = Tw/(1/2)psus? 
specific heat at constant pressure and constant 
volume, (B.t.u./Ib. °R.) 
thermal conductance factor, [Nus/(Prs)'/* 
(Res)' "line / [Nuts (Pr5)'! 3 Re;)' 4liM =o 
kinetic theory function relating thermal con- 
ductivity, viscosity, and specific heat at con- 
stant volume 
acceleration due to gravity, 32.2 ft. per sec. per 
sec. 
G(“), g(n) shearing stress parameters, r = G(u) V/ x, 
G(u) = V ppurl Vhe)/2-¢, g(t) = dn/dt 
wand g, = dg/dt 


or g(t) 
Ln» Rt <i dg dn _ 
local heat-transfer coefficient, he = [R(OT/Oy) |w+ 
(T; — Tw), (B.t.u./hr. ft.2 °R.) 
specific enthalpy, (B.t.u./Ib.) 
molar enthalpy, (B.t.u./mole) 
molar heat of formation, (B.t.u 
mechanical equivalent of heat, 778.2 ft.lbs. per 
B.t.u. 
enthalpy parameter, 7 = h/he 
jn = dj/dyn and j, = dj/dt 


mole) 


pressure equilibrium constant |fi.e., (Kp)o, = 


po?/ por] 
thermal conductivity, (B.t.u./hr. ft.2 °R./ft.) 
Boltzmann constant, kg = 1.381 X 107!* ¢.g.s 
mass of molecule, m = W/(6.02 X* 107%) 
concentration, (Gm. per cu.cm.) 
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p = pressure 

p = partial pressure 

prt = total pressure 

QO = heat transmission per unit area, (B.t.u./hr. ft.?) 

q = dynamic pressure, g = (1/2)psu;5*, (Ibs. per sq.ft.) 

R = gas constant, 1.986 (cal./mole °K.) or 53.33 
(ft.lbs./Ib. °R.) 

S = equivalent cross section for a collision, Sy: = 
1012", etc. 

r = absolute temperature (°R.) 

Ya = plate temperature 

Ts = temperature just outside of the boundary layer 

Tp = reference or base temperature 

fi = insulated plate temperature without radiation 

7. = constant in Sutherland formula (=209°R.) 

t = independent variable of integration, y = 

1 

V QupxX/ppV he f, (u/pp) dt 

u,v = velocity components in the boundary layer 


parallel to the x and y directions, respectively, 


(ft. per sec.) 


U5 = velocity outside boundary layer, (ft. per sec.) 

W = molecular weight 

w = dg dyn 

x,y = Cartesian coordinates; x = distance along the 
plate from the leading edge, (ft.); » = dis 
tance normal to the plate, (ft.) 

= space variable, z = y/+/x 

a, B = number of moles of Ne, Oo required to produce 
1 mole of dissociated products 

7 = ratioof specific heats, y = Cp/¢r 

é = weighting factor for combining viscosities of gas 

mixtures, fi = (Si2e/S))V [1 + (1/me)]/2, 

fo = (Sa /S2)V [1 + (me/m)|/2 

n = velocity parameter, 7 = u/V//p 

bu = absolute viscosity, (Ib.-sec. /ft.*) 

p = density, (lb.-sec.?/ft.*) 

a = molecular diameter, (cm. ) 

T = local shearing stress, r = u(Ou/Oy) = 
V ppuplhe) /*/2x-¢, (Ibs. per sq.ft 

viscosity exponent of temperature, uw «= 7“ 
ff) = constant in Sutherland formula ( =0.522) 


Dimensionless Parameters 


V = Mach Number, (u/a) 


Vu = Nusselt Number, (/.x/k) 
P? = Prandtl Number, (uc,/k) 
Re = Reynolds’ Number, (pux/p) 
(pu/pyu,) = density-viscosity function 


Subscripts 


wW = (pw, Ty, etc.) value at the plate 

6 = (ys, 75, etc.) value at the outside of the boundary 
layer 

B = (pu,, 7p, etc.) value at the reference or base tem- 


perature, B 


D = denotes dissociation 


(I) INTRODUCTION 


oe OF THE LAMINAR BOUNDARY-LAYER EQUA 
TIONS for a flat plate in compressible (air) flow, 
taking into account temperature variation due to con- 
duction and frictional heating, have been studied by 
many investigators. However, these analyses always 
assumed a constant value of Prandtl Number and an 
approximate expression for the viscosity variation with 
temperature. Thus, the applicability to a wide Mach 


Number range is questionable. von Karmén and 
Tsien’ treated both the insulated plate and the hey. 
transfer cases with Mach Number as high as 10, They 
assumed Prandtl Number invariant with temperature 
and viscosity exponents of 0.76, where the viscosity 
exponent, w, is defined by uw « 7°. ) 

Brainerd and Emmons? examined the insulated plat 
case with Mach Numbers up to 3.16 and used Py ~ 
0.733 and w = 0.768. Crocco*® obtained solutions fo; 
both the insulated plate and heat-transfer cases, as 
suming Pr = 0.725 and w = 1.5, 1.00, 0.75, and 0.59 
Hantzsche and Wendt? also treated the insulated plat 
and heat-transfer cases with Mach Numbers up to 1() 
They examined two combinations of air properties 
optional Prandtl Number and w = 1.0, and optional ¢ 
and Pr = 1.0. Cope and Hartree’ used Pr = 0,76. 
w = 0.89 in their treatment of the compressible bound 
ary layer. 

With a greater emphasis on hypersonic speeds at the 
present time, solutions of the laminar boundary-layer 
equations are needed with proper air property-temper- 
ature variation to cover the higher Mach Number 
range. The present analysis covers both the insulated 
plate and heat-transfer cases with Mach Number 
varying from 1 to 20. The laminar boundary-layer 
equations are transformed and reduced to forms suit 
able for solutions with analog computing machines and 
for the employment of experimental variation of air 
properties (u, Pr, etc.) with temperature. At tem 
peratures above which experimental data are available 
(1,800° to 3,000°R.), the variation of air properties are 
extended with the aid of kinetic theory. The air is 
assumed to dissociate to equilibrium in the high 
temperature region, altering the variation of the air 
properties with temperature. The solutions are ob 
tained with a differential analyzer. 

Since the completion of the present work, solutions 0! 
the laminar boundary-layer equations have been made 
by at least three other investigators covering the wide 
Mach Number range. Using Crocco’s method of hand 
calculation, Van Driest'! obtained results for Pr = 0.75 
and a Sutherland law of viscosity-temperature vari 
ation, @ = 0.505, for the heat-transfer plate and the 
insulated plate. A constant value of specific heat at 
constant pressure was used. Young and Janssen! 
solved the reduced boundary-layer equations in integral 
form of the present analysis with the assumption that 
no dissociation occurred. Using the differential an 
alyzer, experimental values of air properties extended 
by the kinetic theory at higher temperatures were intro- 
duced. The Prandtl Number and_viscosity-density 
variations with enthalpy are given by Figs. 6 and 
for no dissociation. 

The results of the several investigators for the in 
sulated plate skin-friction coefficient are shown in Fig 
10, along with the present results for dissociation 
Further comparisons showing the effect of dissociation 


are made with the Young and Janssen results for the 
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Fic. le. 


skin-friction coefficient with heat transfer at the plate 
(Fig. 11), the heat-transfer coefficient (Fig. 12), and 
heat transmission per unit area (Figs. 14 and 16). 


(II) THe THERMODYNAMIC AND PHYSICAL PROPERTIES 
OF AIR 


To be able to solve the reduced boundary-layer equa- 
tions derived in Section (III), the air properties y, c,, 
k, and h must be known as a function of temperature. 
The advantage of the present method of analysis lies in 
the fact that true variation of air properties with tem- 
perature can be utilized without recourse to approxi- 
mate expressions. This is possible because the numeri- 
cal solutions employing analog computing machines 
can absorb graphical data. In terms of the dimension- 
less parameters, the required graphical data are: (1) 
Prvs.j; (2) pu/paurn vs.j; and (3) 7 vs. 7. 

In the high Mach Number range, the temperature 
in the boundary layer can attain such high values as to 
cause dissociation of the air. The properties of the 
dissociated air differ considerably from those of the 
undissociated air. In the present analysis, the rate of 
dissociation is not included, the boundary-layer char- 
acteristics being obtained for equilibrium dissociation of 
the air. The equilibrium composition of air computed 
according to Part (B) of this Section is shown in Figs. 
la, 1b, lc, 1d, and le for several pressures. The steady- 
state condition where the dissociating products of air 
are in equilibrium at all temperatures constitutes a 
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limiting condition. It is recognized that the actyal 
state of the boundary layer in the hypersonic velocity 
régime may lie somewhere between the undissociateg 
and the equilibrium dissociation condition. Also, jt 
might be argued that there is a certain amount of dif. 
fusion of the dissociated products into cooler adjacent 
layers with accompanying recombination. This would, 
no doubt, have an effect on the Prandtl Number. Hoy. 
ever, again, this would be a rate process requiring modi- 
fication of the energy equation. The results, taking 
these rate processes into account, are expected to lie 
in between the results of the steady-state equilibrium 
condition treated in the present analysis and those for 
no dissociation referenced here. Since the results for 
the two limiting conditions so far as skin friction and 
heat transfer of a plate at moderate temperature is con- 
cerned are not appreciably different [see Section (VI 
the rate processes can be neglected. 

For the discussion of air properties as functions of 
temperature the region of temperature below dissoci- 
ation will be treated first, after which the region of 
dissociation will be considered. 


(A). Properties of Air at Lower Temperatures 


In the temperature range below the point where dis- 
sociation occurs, the heat of formation does not enter 
into the determination of properties of air. The vari- 
ation of enthalpy of air with temperature is determined 
from data on the air components O, and N». The 
molecular weight and partial pressures of the com- 
The plot 
of molecular weight for temperatures below dissociation 


ponents are assumed constant in this range. 


is included in Fig. 2. The specific enthalpy is com- 


puted by the expression 
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h = SHip./Wp, 


where H/, is the molar enthalpy based upon National 
Bureau of Standards data® and spectroscopic data giv- 
ing partition functions from which enthalpy is com- 
The thermodynamic 
p; is the partial pres- 


puted by the method of Bethe.’ 
data are compiled in reference 8. 
sure, p, is the total pressure, and W is the molecular 
h is plotted as a variation with temperature 
in Fig. 3. The dimensionless enthalpy parameter j = 
h hg as a function of JT for temperature below the range 


weight. 


of dissociation is included in Fig. 4. 


(hg = 96.8 B.t.u. per Ib. for Ty = 400°R.) 


The specific heat at constant pressure can be obtained 


as 


Cp = dh/dT 


Dp 


In forming the quantities Pr and pu/ppus, experi- 
mental values are selected when available and the ki- 
netic theory is used as a guide to obtain values for higher 
Pr is known experimentally for tem- 
Above this temper- 


temperatures. 
peratures as high as 2,000°R.° 
ature, Pr is obtained by applying the relation" 

f = k/pe, . 
to the components of air, where f = 1.97 for Nz and 
f = 1.92 for Op. 

Pr, expressed as 


Pr = ue,/k = v/f 


is computed for each component using the enthalpy 
data for the determination of y. Pr for undissociated 
air is estimated at the higher temperatures using these 
The variation of Pr with T is plotted in Fig 
5, and its variation with j is included for temperatures 
The experimental vis- 


results. 


below dissociation in Fig. 6. 
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cosity variation of Vasilesco'! is used to a temper- 
ature of 3,240°R. From 3,240°R. to the value where 
dissociation occurs, the relation 


M/ Mer =3,24°R) = (7 3,240)°° 


is used. The variation of absolute viscosity of air for 
the temperature range below dissociation is included in 
Fig. 7. The parameter, pu/pgus, for air below the dis- 
sociation temperature is included as a variation with j 
in Fig. 8, the density variation being a reciprocal func- 
tion of 7 in accordance with the perfect gas law. 

Values of the quantities h, 7, Pr, u, and pu/prus as- 
suming no dissociation are shown for reference in Figs. 
3 through 8 in the higher temperature range where dis- 
sociation can occur |see Section (VI) and also refer- 
ence 15]. The viscosity variation with temperature 
for the undissociated air at the higher temperatures 
referenced in Fig. 7 is calculated according to the 
Sutherland formula, 


Be r) _». 366 
Mp s VT/Tz + (6/0 T/Ts) 
where 6 = 7;/T, = 0.522. 


(B). Properties of Air at Higher Temperatures 


At the higher temperature where dissociation occurs, 
the chemical composition of the air had to be deter- 
mined in order to evaluate the air properties by the 
kinetic theory. The chemical process of dissociation 
essentially is 0.79p,N2 + 0.21p,O2 + heat > piN2 + 
pxO2 + psN + psO + psNO where p; are the partial 


pressures of the various constituents. For dissoci- 


ation equilibrium, the following five equations are 
solved simultaneously : 


Chemical balance: 


2p1 + ps + Ps me 
0.79), 


2po + pa + Ps 
0.21p, 


Dalton’s Law of partial pressures: 
Zpi/pr = | 
Equilibrium equations: 


(K,)o, _ 


Po’ Po, 
(Ky), me pr? Py. 
(K,)wo = pno (pw,)' *(Po.) ; 


where the A, are the pressure equilibrium constants’ 
and p, is the total pressure. Partial pressure of the air 
components were computed for several total pressures 
(10-4 to 1 atm.) covering a temperature range up to 
14,000°R. The results are presented in Fig. 1. The 
molecular weight of the mixture of the dissociated prod- 
ucts, 2W;p;/p,, is shown in Fig. 2, W; being the molec- 
ular weight of the 7th component. 

The specific enthalpy including the heats of forma- 
tion, shown in Fig. 3 for pressures 0.0001 to | atm., is 


calculated by 
hp = (2Hip; + 2H/pi)/=Wip, 


where //; are the molar enthalpies and //,/ are the heats 
of formation based upon National Bureau of Standards 
data® and spectroscopic data giving the partition func- 
tions from which enthalpies are calculated by the 
method of Bethe.’ 

These data are compiled in reference 8. 
sionless enthalpy parameter j = (h/hg)» for dissociated 
air as a function of 7 is presented in Fig. 4. The sub- 
script D used here denotes dissociation. The specific 
heat at constant pressure in the dissociation range is 


The dimen- 


expressed as 


Cp, = (dh/dT)p 


’D 


and the parameter Prp for dissociated air is computed 


by 


—_ > - - 
Prp = Pr(cy,,/Cp) 


where Pr is the Prandtl Number for undissociated air 
(Fig. 5). An assumption herein implied is that the 
ratio u/k is relatively unaffected by dissociation, since 
both viscosity and thermal conductivity are transport 
properties. Inasmuch as the effect of pressure is small, 
Pry vs. j (Fig. 6) for p = 0.01 atm. is used for all pres- 
sures. 

The plot of (pu/pzus)p VS. j Shown in Fig. 8 for equilib- 
rium dissociation of air is also essentially independent of 


pressure. The density ratio is determined by 


p/pp = (W/Wsa)-Tp/T 
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LAMINAR BOUNDAR 


In obtaining the viscosity ratio for dissociated air by 
kinetic theory at the higher temperatures, the method 
is as follows: The viscosity expression for the physical 


model that assumes uniform elastic spherical molecules 


of diameter o 1S 


— 5 { tant) ; = 0.1792 Com!) 
160° a” 


T 


where kg is the Boltzmann constant and m is the mass 
of molecule. Treating all the monatomic species as a 
group, the viscosity is computed according to this 
Similarly, the viscosity is evaluated for the 


model. 
The resultant viscosity of the mix- 


diatomic species. 
ture of monatomic and diatomic constituents is hence- 


forth obtained by the relation’ 


eal Me 
i 1 + §12(2/m1) L + &1(1/M2) 
where 

; Si [1 + (m/m)]'”* 

gr _ S; V/2 

: Sx [1 + (me2/m)]'”? 

fo "<< /2 

Ny, Nx concentrations 

Mm, M, = average masses of molecules of monatomic 
and diatomic species, respectively 

Si = mo," = cross section of collision of the first 
kind with the first kind (monatomic) 

Se = mo.” = cross section of collision of the 
second kind with the second kind (dia- 
tomic) 

Si = S» = row? = cross section of collision of 
the first kind (monatomic) with the 
second kind (diatomic) 

01, 02 = diameters of first and second kinds of 
molecules 

1 = (1/2) (o; + o2) = mean diameter 


A graph of the viscosity of air for equilibrium dissoci- 
ation is given in Fig. 7. 
SUITABLE FOR 


(III) DerrivaTION OF EQUATIONS 


ANALOG COMPUTING MACHINES 


The laminar boundary-layer equations for compres- 
sible flow over a flat plate, including temperature vari- 


ations, are: 


Conservation of momentum: 


Ou 1 Ou 2 ( o) 
” — oy = oy ‘ oy 


Conservation of mass: 


[O(pu)/Ox] + [O(pv)/Oy] = O 
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Fic. 8. Density-viscosity parameter vs. enthalpy parameter 


Conservation of energy: 


oh oh oO ~*) on 
= Ox ie oy 7” mG oy bie & 

Zero pressure gradient along the plate is assumed. The 
enthalpy (or heat content), 4, includes the heat of for- 
mation of dissociated products at the higher temper- 
atures where dissociation occurs. 

Since these equations are partial nonlinear differ- 
ential equations, their solution is most formidable. To 
obviate the difficulty, an assumption is borrowed from 


incompressible boundary-layer theory. This assump- 


tion is that u = f(y/+/x)—i.e., the velocity—and, 
therefore, 7, h, u, and k are constant along parabolas 
through the origin. Thus, 

u= u(z),v = v(z), 7 = T(z), wu = w(s),... 


where z = y/+/x. The three partial differential equa- 
tions can be shown to reduce to three total differential 
equations in z alone. This results from the 
that the equations obtained in transforming from 
independent variables x, y to x, z are of such form that 
the terms containing z and its derivatives do not con- 
With the velocity and air properties taken as 
Upon 


fact 


tain x. 
functions of z only, the x terms can be dropped. 
transforming the independent variable from z to the 
velocity variable u and introducing a function G(u) = 
u(du/dz), the combination of momentum and energy 


a ’ 
( +G 


7 
uc, du 


equations yields 


dG dh d ( k 


dudu = du 
Combining the momentum and mass equations yields 


G(d°G/du?) = —(1/2)pyu 


Up to this point the derivation is identical with that 
of Hantzsche and Wendt.‘ These two equations are 
now made dimensionless, with the following dimension- 


less variables: 
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Velocity: 
"n= u/© he 


Shearing stress: 
g(n) = G(u) Vv (1 2) pps hip) 
Enthalpy: 
j=h/hp 


where the subscript B indicates the value of the air 
properties at a base or reference temperature of 400°R. 
The two equations then become 


dg dj 1/{gdj 
gdj_ ¢ (: <2) 4 
dndn  dn\Prdn 


g(d*g/dn*) = —(ppu/paun)n 


where Pr = uc,/k. 

Pr and (pu/prun) contain all the physical properties 
of air including the effects of dissociation. These quan- 
tities can best be expressed as functions of 7, a thermo- 
dynamic property of air also affected by dissociation 
[Section (II) }. 

The above equations must be transformed to a new 
independent variable to eliminate an infinity value at 
the upper limit (i.e., |/g — ©), since this condition is 
machine methods. The 


unsuitable for computing 


transformation eliminating this singularity is 
g = dn/dt 
It is also convenient to define a function w by 
dw/dt = (pu/pspes)n 


The two equations then become 


_ dG dfil a) 4 92 
dt dt\Prdt) | * 


d*j p dj ‘ (1) Pr) 4 | 
= — ri W z2 
dt? dt‘ \dt) dt 6 


dg/dn = —w 


or 


and 


since w = 0) when y = 0. 
The five first-order equations to be solved simul- 
taneously by the differential analyzer'’ are, in integral 


form, 
n= Sgdt 
w= S (pu/peun)n dt 
g = —fwdr 


lj ‘dj [1 
—S Prd fwdj + | nal ) + fg in| 
J Gt NIT 


S (dj/dt) dt 
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The equation relating the original independent vay. 
ables, x and y, to the independent variable of integra. 


9 is I/» ‘/ 

2upx Me 
la ( = ) | dt 

PBV hp 9 MB 


The five boundary conditions to accompany the jn. 


tion, /, is 


tegral equations are all considered at ¢ = 0, correspond- 
ing to y = 0-—1.e., at the plate surface. They are: 

(1) » = 0; the velocity components u and 7 are zero 
at the plate. 

(2) w= —dg/dn; the gradient of the shearing stress 
is zero at the plate. 

(3) J = Jw, the enthalpy of the air next to the plate 
determined by the wall temperature, which is to be 
fixed for each set of calculations. 

(4) g = gy, the shearing stress at the wall fixed 
for each set of calculations. 

(5) dj/dt = [dj/dt],, the enthalpy gradient deter- 
mined by the temperature gradient or heat-transfer 
condition at the wall. This value also is to be fixed for 
each set of calculations. For an insulated plate with 
no radiant heat transfer, [dj/dt|, = 0. 


(IV) COMPUTING MACHINE METHODS AND SOLUTIONS 


The five first-order equations in integral form [Sec- 
tion (III)] are solved by assigning the boundary con- 
ditions at the plate surface ¢ = 0 and allowing the ma 
chine to proceed with the solution as / goes to ~. The 
terminal values of temperature at the outside of the 
boundary layer and local Mach Number are thereby 
determined by the solution. It would, of course, have 
been desirable to have been able to specify nominal 
values of 7; and J/;, as well as 7), and temperature 
gradient at the plate, in order to standardize the re- 
sults. However, since interest was primarily in skin 
friction and heat transfer and not in the distribution of 
shear, temperature, and velocity through the boundary 
layer, an alternate method was available. Values of 
skin-friction and heat-transfer parameters can be ob- 
tained for nominal values of the conditions at the plate 
and the outside of the boundary layer by making a 
series of runs, the results of which are plotted. Hence, 
for a chosen temperature outside the boundary layer, a 
set of runs was made whereby four of the five initial 
boundary conditions were held constant and the fifth 
one varied. Initial parameter estimates were made 
which gave values of j; covering the range from 0.5 to 
4, and a cut was taken at j; = 1(7; = 400°R.).T 

For the insulated plate solution, (dj/df),, = 0, the 
shearing stress at the plate, g,,, was made the variable 
Plotting g,. vs. j; and 7 
For the 


parameter for a set of runs. 
. >. . - . 
VS. Js, Zw and ns; were obtained for j; = 1. 
t All calculations of friction and heat-transfer coefficients 10 
this analysis are based upon 73 = 400°R. However, the data 
can be used to obtain results for other values of 75 by taking other 
cuts in the range 200° to 1,600°R. 
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TABLE | 
ind 7; for 76 = 1 for the Insulated Plate, 
(dj/dt)w = 0 


Values of Jw, Zw, 


Lu 76 
3.0 1.37 2.23 
6.0 2.49 3.57 
12.0 3.96 5.15 
25.0 4.98 6.36 
50.0 6.71 8.06 
75.0 7.63 8.98 
125.0 9.71 10.98 
137.5 11.038 12.18 
145.0 12.17 13.33 
TABLE 2 
Values of jw, (dj/dt)w, Zw, and ns for js = 1 for the case of Heat 


Transfer at the Plate 


g de (dj/dt)u 5 
0.30 | 0.75 0.74 
0.69 l 0.34 1.30 
2.12 | 2.33 2.82 
2 1.59 2.93 
3 1.04 3.02 
12 —5.5 3.36 
25 —25.2 3.44} REACT 
50 —60.8 3.55 
3.5 l 5.55 $1] 
3 3.70 1.32 
6 1.82 4.48» REAC 
25 —18.7 4.80 
50 —47.7 5. 17 
$4.95 3 8.00 5.60 
6 6.01 5.79 
12 4.18 6.04 
1S 3.80 6.16 
50 —39.5 6.32 " 
75 —104 6.73f REAC 
6.5 l 16.8 6.60 
} 13.7 6.85 
6 10.9 7.00 “ee 
12 12.5 7.39/ REAC 
25 15.2 7.45 
795 — 58.4 1.43 
8.20 ; 21.8 8.24 
6 19.2 8.41 
12 24.1 8.82 
5 47.5 9.12 
1) 34.7 9.37 
75 14.5 9.43; REAC 
10.0 l 36.7 9.34 
6 28.0 9.70 
12 40.0 10.12 
18 75.3 10.24 
25 80.3 10.33 suee 
34 54.8 10.37; REAC 
410) 62.9 10.45 
50 wa.8 10.55 
60 76.5 10.59 
SO 62.3 10.69 
100 2 10.75 ° 
11.50 3 10.5 10.71 
6 37.5 11.09 
12 §1.2 11.55 
25 119.0 11.81 
50 131.2 12.08 
75 145.8 12.17 
100 109.3 12.39 
14.50 ; 61.6 13.28 
6 57.7 13.49 
25 187.2 13.41} REAC 


t Additional solutions obtained for the equilibrium dissociation 
case on the Reeves Electronic Analogue Computor are designated 
‘REAC.”” Experience showed the REAC to be less accurate 
but much more rapid in running solutions. Its results used as 
eeeventary data to the differential analyzer results are satis 
actory 
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case of heat (dj/dt)» # O, (dj/dt), was 


made the variable parameter, keeping the plate tem- 


transfer, 


perature and the shearing stress of the plate fixed at a 
constant value. (dj/dt)» vs. j; and 75 vs. js were plotted 
to obtain values corresponding to j; = 1. 

Values of the dimensionless parameters for j; = | 
are shown for the insulated plate and the plate with 
heat transfer in Tables 1 and 2. Values of 7y, Tw, 
(OT/Ov)~, ws, and 7; can be obtained from the dimen- 
sionless variables. 

Typical graphical results of boundary-layer distribu- 
tions obtained by using the differential analyzer are 
shown in Fig. 9. Both the insulated plate case and that 
of heat transfer to the plate are represented, as well as 
low and high Mach Numbers. Velocity and temper- 
ature distribution u vs. y//x and 7 vs. y VV xf can 
be obtained from the graphical results, although the re- 
sults will not be for nominal plate and outside boundary- 
layer conditions. The transformation from the vari- 
able ¢ to the variable y, which is the distance normal 
to the plate, is made according to 


9 .\} 
~MBX LM 
( - ) | dt 
ppv hp 0 MB 


Prior to solving the boundary-layer equations, the 
This was 
In the 


¥y = 


accuracy of the machine setup was checked. 
accomplished by two simple comparison tests. 


first test, made at J, = 2 and M, = 5, the relation 
Cys V Res = V2 Sw/(ns) ‘> = 0.664 for the insulated 
plate for Pr = 1 and w = 1| was duplicated. Another 


check, made at a high Mach Number value with vari- 
able air properties, was shown to agree with hand- 
computed results 
(V) RersuLts—SKIN FRICTION AND HEAT TRANSFER 
The differential analyzer solutions are used to obtain 
the skin-friction and heat-transfer parameters in a form 
+ Since the velocity and temperature distributions can be ob- 
tained for each set of calculations, the boundary-layer thickness 
6, displacement thickness 6*, and the momentum thickness 6** 


can also be determined 
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applicable to the field of aerodynamics. The local fric- 


tion coefficient is treated first. 


(A) Local Friction Coefficient 


Defined as the ratio of the unit shearing stress at the 
plate to the dynamic pressure, the local skin-friction 
coefficient is expressed as 


ne V Re, = Tw pg sX 
(1/2)psus"L Ms 
Since 
To = V paun(hs)’ */ 2X Lu 
The local skin-friction coefficient is related to the 


mathematical parameters as follows: 


; ie 3/s V/s 
Cys V Res = V2 Bw/ 05 °° (05Hs/ Paes) 


where air properties are based upon the condition at 
the outside of the boundary layer as designated by the 
subscript 6. 

The variation of c,;; ~/ Re; with M, is computed from 
the results of this analysis using Tables 1 and 2 for an 
ambient temperature 7; = 400°R. The results for the 
insulated plate and the plate with heat transfer with 
the plate enthalpy maintained at several constant values 
are shown in Figs. 10 and 11, respectively. See Section 
(VI) for comparisons with the results of other investi- 


gators. 
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Mach Number at the edge of the boundary layer is 
determined by 


u V 32.27h 
M;, => , = : Bins = L.5OS5S8 1 
@ V32.2yRTRV j; Vi; 
where 

J = W782 tt lbs. /3B-t:u. 

he = 96.8B.t.u./Ib. 

R = 52.3 ft.lbs./tb. °R. 

Y = C)/cC, = 1.405 

Tn = 400°R. 


(B) Heat-Transfer Coefficient 


The expression for convective heat transfer through 
the boundary layer can be written as 


oT —— 
Q = ke(.) = h(i — Te) 
Vw 


which defines the local heat-transfer coefficient in terms 
of the insulated plate temperature, 7%. 


Q = heat transmission, (B.t.u./hr.ft.’) 

k, = thermal conductivity, (B.t.u./hr.ft.? °R. /ft 
h. = heat-transfer coefficient, (B.t.u./hr.ft.2 °R. 
7, = plate temperature, (°R.) 


The temperature gradient at the plate in terms of 
the mathematical parameters is expressed as 


(OT /Oy)» = (dT /dn)w (dn/du) (Ou/Oy), 
where 


(dT/dn)w = (Ap/CpwSw) (dj/dt) x 





dn/du = 1/Vhp 


(o") Tx gw lppuBlp 

oy u Ky 

When these quantities are introduced into the heat- 
transfer equation to form a parameter, Nu;/¥V Pr 


Mic 2x 


V/ Re;, the following expression results: 


Nu; l (Pr;)’ (@) Cpr l x ' 
\ Prsv/ Re; (T; _ re (Pri) dt u ¥ Crs V 2 
] seh 
Vv Nal. Pos 


which is proportional to the heat-transfer coefficient 
For incompressible 


for a constant value of velocity. 
laminar flow over a flat plate with no dissociation, the 


expression is 
Nu;/V PrsV Res| we = 0.331 


which also holds for compressible flow when local values 
of the properties and variables are taken. To show the 
effect of dissociation on the heat-transfer coefficient 


the hypersonic flow régime, the ratio 
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is computed for several plate enthalpy values using the 
data of Table 2. The variation of this thermal con- 
ductance factor with V7; is shown in Fig. 12 for several 
values of j, (or 7.) and p. The laminar heat-transfer 
coefficient, h., for dissociation 
by applying the above factor as a correction to h, 


can be computed 
for laminar incompressible flow with no dissocia- 
tion. 

Insulated plate temperature, 7), is obtained from the 
results of Table 1 giving values of 7,, and y;. These are 
reduced to 7; and u;, respectively, using Fig. 4 and the 


relation 
us = V 32.2Jhe ns = 1,547n, 


The results are plotted as a temperature rise above 
T;, AT; vs. us in Fig. 13. Shown for comparison pur- 
poses are plots of AT; vs. u,; for no dissociation and for 
dissociation equilibrium the insulated plate 
temperatures are computed from the energy relation 
The latter result for 


where 


with a recovery factor included. 
dissociation equilibrium takes into account the heats 
of formation of dissociated products but does not take 
into account the variation of temperature through the 
boundary layer. It would be true if the air were stag- 
nated en masse. 
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(C) Heat Transmission per Unit Area 
The expression, Q = h,(7; — T,,), giving the convec- 


tive heat transfer through the boundary layer, is ap- 
plicable to flow in the hypersonic flow régime provided 
h, and 7°; are modified to take into account dissociation. 
A more direct method is to compute QO simply by ex- 
I ply )) 
pressing k, (O7'/Oy), in terms of the mathematical 


parameters of this analysis. 
the relation is 


QO = k,(0T/dy), 


| 
Ppsur hgGJ (dj dt 


2x PY, 


= 3.6 )OGhz 


(dj dt), 


= 2.972 X 10° V p/x 
Vy 


Using engineering units, 


(B.t.u./hr.ft.”), where p is given in atmospheres and x 


in feet. The variation of (dj/dt), with 


the velocity 


parameter 7; is plotted for both moderate and high 
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Fic. 13. 


values of enthalpy at the plate in Fig. 14. The value for 


Pr, is taken from Fig. 6. 


(VI) CONCLUDING REMARKS 


It is interesting to compare the skin-friction and 
heat-transfer results obtained in this analysis which 
assumed dissociation equilibrium of the boundary layer 
air with the results of others for which no dissociation of 
air is assumed. In each case the boundary layer on the 
flat plate is laminar and the flow is compressible. 

The skin-friction results for the insulated plate by 
several investigators with various assumptions of air 
properties noted are shown in comparison with the re- 
sults obtained here in Fig. 10. Of particular interest 
are the results of Young and Janssen” for Mach Num- 
bers approximately | to 20, wherein true air property 
variation with temperature without dissociation was 
assumed. Their results agree well with the results of 
this study for the skin-friction coefficient on the insu- 
lated plate. Agreement is also good between the disso- 
ciation and no-dissociation results for the skin-friction 
coefficient of the plate with heat transfer occurring for 
values of enthalpy at the plate j, = 3 and 
The comparison can be 


moderate 
jw = 6 as shown in Fig. 11. 


TIL at 
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oZ 


made equally well by plotting the shearing stress param 
eter g,, VS. ns resulting in one curve for both dissociatio; 
and (Fig. It can therefore be 
stated that dissociation has a negligible effect on the 
local skin-friction coefficient for the insulated plate 


no dissociation l5a). 


and the plate with heat transfer occurring for plate tem 
peratures 400 to 2,500°R. 
between the dissociation and no-dissociation curyes 


However, some divergeng 


occurs for higher values of enthalpy at the plate (j,, = 95 
and j. = 75). 15b. The effect 
of dissociation is to increase the friction coefficient jy 


This is shown in Fig. 


this case. 
The effect of 
efficient is observed by comparing the thermal con 


dissociation on the heat-transfer oq 


ductance factor, F, computed for no dissociation but 


variable air properties’ with F computed from this 


analysis. This is done in Fig. 12, where a few points 
are plotted to show that the ratio / remains close to 
unity through the Mach Number range for both moder. 
ate and high values of enthalpy at the plate for the no- 
dissociation assumption. The in 
of heat-transfer coefficient with Mach Number for th 


increase the value 
dissociation solutions thus appears to be caused entirely 
by the dissociation process and not by the variability 
in the properties of air, u, k, and c,, etc., when local 


Dissociation causes a_ significant 


T; (Fig. 


values are taken. 
lowering of the insulated plate temperature, 
13). 

on heat 


To observe the effect of dissociation the 


transmission, Q, it is necessary only to compare values 


of (dj/dt), for the two conditions. For moderate 
values of enthalpy at the plate, 7, = 3 and j,. = 6, one 
curve is common in plotting (dj/dt)» vs. ns for dissoci 
ation or no-dissociation results (Fig. 14). Hence, the 
heat transmission, Q, is equal for the two cases. How 
ever, for larger values of j,, (1.€., jy = 25 and jy, = 75), 
the agreement is no longer evident (Fig. 16). At the 


extremely high plate temperatures, the effect of dis 
sociation on Q is one of cooling at lower Mach Number 
and heating at higher J/. 
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Fic. 14. Gradient of enthalpy parameter at plate as a func 
tion of velocity parameter (A) with and without dissociation 
(jw = 1,jw = 3, and j» = 6) and(B) with'dissociation (jw = 12, 
jw = 25, and jw = 50). 
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TABLF 3 
25 

g (dj/dt)u 5 
0.30 —6.40 0.84 
0.94 —8.63 2.00 
1.46 —9 30 ° 70 
1.85 —9.70 3.25 
2.58 —9.15 1.12 
3.18 —8.90 1.71 
4.83 —5.75 6.31 
7.438 3.45 8.65 
8.90 9.63 9 88 
10.21 16.50 10.93 
11.42 23.50 11.82 
13.14 35.20 13.19 
15.00 $7.00 14.61 
= 75 

g dj dtu "5 
0.30 —16.78 lo 
0.94 — 23.19 2.3 
1.46 — 26.00 3.04 
1.85 —27.50 3.58 
2.58 —29 84 1.55 
3.18 —30.35 5.24 
4.83 — 30.45 7.07 
5.80 — 28.85 7.95 
7.43 — 25.25 9.59 
8.90 20.30 10.90 
10.21 —14.85 11.96 
11.42 —9 30 12.93 


The criterion for determining if dissociation is to 
affect the skin-friction and heat-transmission param- 
eters is as follows: The temperature at the plate must 
be high enough for dissociation to occur there, conse- 
quently affecting the shearing stress and gradient of 
enthalpy or temperature at the plate. cy, for the in- 
sulated plate appears to be an exception to this rule, 
since it is unaffected by the dissociation occurring at the 


plate. 


APPENDIX (A) 


Laminar boundary-layer characteristics for no disso- 
ciation but variable properties of the air were obtained 
for high plate enthalpy values from later solutions? of 
the boundary-layer equations. Since these results have 
not been reported, tabulations of (dj/df),, go, and 75 
are given for j, = 25 and j, = 75 for js = 1 1n Table 3. 
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ABSTRACT 


Whirling of a bladed disc mounted on a cantilever shaft is dis- 
assuming that the elastic blades have a highly idealized 
form. Results show that whirling cannot take place, except 
under restricted conditions, if shaft possesses symmetrical elas- 
When whirling can occur, blade flexibility has 


cussec 1 ’ 


tic properties 
in appreciable effect on some of the whirling velocities. 


NOTATION 


center of gravity of disc 


c.g 

Ui = mass of rigid or bladed disc 

WA = mass of disc in bladed disc 

M = total mass of blade tips 

] = diametral moment of inertia of rigid disc or bladed 
disc (static) 

Ly = length of cantilever shaft 

L = radius of bladed disc 

| = length of blade 

p = radius of circle described by disc c.g 

B,, By = shaft flexibility 

wp = angular static natural frequency of blade 

w = angular velocity of disc 

¢ = whirling angular velocity 

8 = ww; —-¢ 

x = inclination of disc axis to OZ 

€ = blade angle 

i = VV -1 

= 2.7182818284 

T, ¥ = angles 


real number 


(1) INTRODUCTION 


eee CONSIDERABLE ATTENTION has been given 


in the past to the problem of whirling of a rigid 
disc mounted on an elastic shaft, little effort has been 
directed to the study of the whirling of a similarly 
mounted rigid disc to whose rim is attached a large 
number of evenly spaced identical elastic blades, as in 
a gas-turbine wheel. The reason for this division of 
attention appears to be due to the fact that, in most in- 
stances in steam turbine and stationary gas-turbine prac- 
tice, the ratios of the lowest static natural frequencies 
of the blades to the static point-mass natural frequency 
of the disc-shaft system are so high that the blades can 
be treated as rigid insofar as whirling is concerned. 
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Whirling of a Bladed Disc’ 


JOHN L. BOGDANOFFT 


Purdue Unwersity 


With the advent of the aircraft gas turbine and axial- 
flow compressor, in which the frequency ratios just 
mentioned are much lower, the problem merits rein- 
vestigation with a view to determining the influence of 
flexible blades on the whirling natural frequencies (or 
whirling velocities). This paper is devoted to such an 
investigation for a bladed disc mounted on a cantilever 
shaft, the blades having a highly idealized form. 

The blades consist of identical massless cantilever 
beams carrying concentrated tip masses. We assume 
that the cantilever consists of a sufficiently wide strip 
of thin elastic material so that the tip mass can be as- 
sumed to move only in the direction perpendicular to 
the plane of the unbent strip. Each blade is mounted 
on the central dise in such a way that the unbent axis 
of a blade lies in the plane of the central disc and is 
coincident with a radial line through the mass center 
of that disc. The blades are evenly spaced around the 
disc circumference, and their number is sufficiently large 
so that the discrete distribution can be replaced with 
sufficient accuracy by a continuous one. 

While it must be admitted that a blade of the type 
just described is not encountered in practice, never- 
theless, an analysis based upon this hypothetical blade 
is permissible at the present time since qualitative 
rather than quantitative information is desired. Hence, 
this paper must be considered as a preliminary report. 
The effects of such blade factors as continuous mass 
distribution, twist, noncoincidence of mass, and elastic 
axis, which are of unknown importance, await further 
investigation. 

The ratio of chief interest is, as has already been 
noted, the ratio of the static blade natural freyuency 
to the natural frequency of a particle on the end of the 
cantilever shaft, the mass of the particle being equal 
to the mass of the bladed disc. Results are presented 
for three values of this ratio—namely, 1, 2, and 4. 
For the system considered, the influence of blade flexi- 
bility on whirling was small when this ratio exceeded 

The remaining ratios encountered in the whirling 
velocity equations are taken to have the following 


values: 


M/M, = 1/4, L/L, = 2 
L/l = 3, 3ko?/2L,? = 1 


Other values could have equally well been chosen. 
However, it is felt that the above values are fairly 
representative of what might be encountered in prac- 
tice, and they also simplified the numerical work toa 
certain extent. 
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velocity denoted by ¢. If w; represents the constant 1! 
p' angular velocity of the disc about its axis, the magni- : 
tudes of the force F and moment N which must be 
exerted at O’ are, respectively, 





Mope’, Io (2w3 = ¢) ex 














D Eq 
O' where the force acts ina direction parallel to O"0 ang Joc 
' the moment is perpendicular to the plane 00/0". Sing is { 
+ the bent shaft exe rts a force and moment at O’ parallel Ws | 
i to the two just mentioned, having magnitudes 
(12B, L,*)p — (6B, L1")x 
(6B, L1")p = (4B, L1)x 
Th 
SS the equations governing whirling motion are fre 
Myp¢? = (12B,/L1*)p — (6B,/L,")x ) 
O I)(2w; — ¢)¢x = (6Bi/Li")p — (4Bi/Li)xf - 
977777 These equations will have a nontrivial solution if, and 
Fic. 1. Equilibrium position of rigid disc on cantilever. only if, 
ma 





The whirling of a rigid disc on a cantilever has been 
included because of its usefulness in comparison pur- 1 
poses and in order to simplify the work that follows. P 
Since this is an old problem, which has been discussed 
by many writers, the treatment is brief and no detailed 
comments are made. The Appendix [with equations 
numbered (A-1), (A-2), etc.] contains a derivation of 
the equations that determine the whirling velocities 
and mode shape for our bladed disc system under the 
assumption that the shaft possesses unsymmetrical 
elastic properties. Little additional effort was re- 
quired to analyze the unsymmetrical case, and its in- 
clusion makes possible the observation that whirling 
cannot take place, except under special circumstances, 
unless the shaft is unsymmetrical. The remainder of 
this paper is devoted to a discussion of the whirling of 
the bladed disc system. 





(2) Ricip Disc ON CANTILEVER SHAFT 


Fig. 1 shows the equilibrium position of a rotating 
disc D mounted on a cantilever shaft S. We assume 
that, when the system is at rest, the axis of the shaft, 
the disc c.g. O’, and the disc axis O’P’ lie on the same 
straight line; S is taken to have a circular cross section 
of constant radius; and gravity is neglected. An 
instantaneous position of the system in motion is shown 
in Fig. 2. The XYZ axes, with origin O at the built-in 
end of the cantilever, are fixed in space, OZ being the 
equilibrium position of the shaft and dise axis; O’Z’ 
is parallel to OZ. The projections of O’ and P’ 
on the OXY plane are denoted by O” and P”, respec- " 
tively. xX fe 

During whirling, O’P’ and S remain in a common \ 
plane through OZ, with p and x equal odanesnere Fic. 2. Instantaneous position of rigid disc on cantilever. Th 
and this plane rotates about OZ with constant angular ‘ angle XOP” should be labeled ¢. 
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- + [)(20:, — ¢'¢ 
L;’ Ly 
Eq. (2) determines the possible whirling angular ve- 
locities ¢ for given w;. The ratio Lx /p, or mode shape, 
js found from either of Eqs. (1) using values of ¢ and 
w, consistent with Eq. (2): 
lix (12B, L*) —_ Moe? 


p 6B,/L;" 


ww 


The introduction of the angular point-mass natural 
frequency of the system, 
a = V3Bi/ ML) 
and the substitutions 
ko? = Io/My, XX = G/a 


make it possible to write Eqs, (2) and (3) in the forms 


3k? W3 
=a (9 ~ —I a ~ 2 fa) ~6 <0 
2L;" a (4) 


Lix/p = (4 — X*)/2 | 
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> "Thi 


Fic. 3b. Mode of motion of rigid disc system. 


Eqs. (4) can also be derived from Eqs. (A-7) and 
(A-S) by putting 6 = Oand letting wz, > @. 

The simplest method known to the author of finding 
the permissible \ for given w;/a is to interchange the 
roles of \ and w;/a—i.e., find the permissible w;/a for 
given \. This is readily accomplished, since the first 
of Eqs. (4) is linear in w;/a. Fig. 3a presents this infor- 
mation for the particular case 3k)? = 2L,°. The mode 
shape, as a function of \, is found from the second of 
Eqs. (4), and the results are given in Fig. 3b. 

It is seen from Fig. 3a that for each disc w; there are 
four distinct angular whirling velocities, two in the same 
sense as w; and two in the opposite sense. These whirl 
ing velocities, when divided by 27, give the whirling 
natural frequencies of the system. 

For detailed comments on the points of interest in 
Figs. 3a and 3b, see a previous paper by the writer.' 


(3) BLADED Disc ON CANTILEVER SHAFT WITH b = O 


Eqs. (A-7) and (A-8) reduce in the present case to 


in(m — 1)oz; = O 

o102 — 6 — (3/32) (1 — 2n)?a;"7 = O 

oi|sin y + (1/8) (1 — 2n)o3 cos y| = 0 (3) 
Lix a  p o1 

p 7 a} ~ 2 cos y — (1/4) (1 — 2n)?o3 sin ¥ 
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Fic. 4a. Whirling velocities of bladed disc system e = 0, 
b = 0, wa/a = 1 


where the ratios have their previously defined values 
and 


o, = 4 — [1 + (sin? €/8) (1 — 2n)@)d? | 
o. = 2 — (1 + (3/4) cos? e@) (1 — 2n)d? (6) 
a3 = (sin 2e/2)Pd? 


We note that if w, > ©, Eqs. (5) reduce to Eqs. (4) 
i.e., the system of this section is dynamically equivalent 
to the one in the previous section if the blades are as- 
sumed to be rigid. The quantity \ is determined from 
the first two of Eqs. (5) as a function of ¢, m, and w;/a;s 
the corresponding y and 1;x/p are found from the third 
and fourth of Eqs. (5), respectively. 
not appear in Eqs. (5), since it has significance only 
when b ¥ 0 

Examination of the first of Eqs. (5) indicates that 
this equation can be satisfied for \ # 0 only if either 
The first of these 


The angle 7 does 


sin 2e = Oor n(n — 1) (1 — 2n) = 0. 
conditions shows that « must equal either 0, 7/2, 7, 
3x/2, etc.; of these values, only the first two represent 
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different equilibrium configurations of the blades with 
From the third of Egs 
(5), it is seen that sin y = 0 for « = 0, 7/2. Thus, y 
equals either 0 or 7; this means that during steady 


respect to the central disc. 


precession the rotating plane through OZ and 0’ also 
contains the central disc axis and the shaft. The 
second condition shows that the first of Eqs. (5) is 
satisfied for all \ and e¢ if m equals either 0, '/, or | 
The case n = 0 is of no technical importance, since it 
concerns the motion of the system when w; = 0; hence, 

With » 
or 1, the second, third, and fourth 


it will be excluded from further discussion. 
equal to either '/2 
equations of Eqs. (5) determine A, y, and Lyx /p as 
functions of «. Since y, in general, does not equal 
zero, S is bent into a space curve during steady motion 
We now shall consider in detail the four cases: (j 


e= 0: (i) « = 90°: (1) 8 = ° and (iv) ” i. 


Case (i): « = 0 


Here, as already mentioned, € equals zero and Eqs, 


(5) become 


010% 6=0, Lix/p = o;/2 7 
where 
o=4-—M ( 
& 
o. = 2 — [1 + (3/4)4] (1 2n)r24 


The heavy solid curves in Fig. 4a represent the dimen 
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Fic. 4b. Mode of motion of bladed disc system: ¢ = 0,5 = 0 
op/a = Ss 
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Fic. 5. Whirling velocities of bladed disc system: € = 0, b = 0, 


wp/a = 2 

sionless whirling angular velocities as a function of the 
dimensionless disc angular velocity w;/a for a disc 
having blades with w, = a. We note that there are 
six possible \ for each w;/a, three in the same sense as 
w; and three in the opposite sense. For ease in discus 
sion, we label these curves Ag, A4, Av, Ai, Az, As, iN going 
from the top to the bottom of the figure. The light 
solid curves determine \ as a function of w;/a when 
w, equals infinity—-i.e., these curves are the same as in 
Fig. 3a. The heavy dashed curves represent the values 
of \ and w;/a@ at which resonant vibration occurs in the 
Fig. 4b gives the ratio L;x/p as a function of X. 

The presence of six whirling velocities for each disc 
velocity, instead of the four encountered in the rigid 
dise system, was to be expected, since we have effectively 
introduced two additional degrees of freedom by per- 
mitting the blades to deflect. 

In those cases examined, \; and Ags are of sufficient 
magnitude so as to be of little technical interest beyond 
hence, 


blades. 


noting that they approach infinity as w, > ©; 
these two whirling velocities will not appear in several 
of the figures pertaining to whirling velocities given 
below. 

We observe from Fig 4a that \; and 2 are not ap- 
preciably different from the comparable rigid disc curves 


Over most of the range of w;/a. Examination of Fig. 


A 
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tb and Eq. (A-9) indicates that on these two curves the 
magnitude of the disturbance that produces motion of 
the blade tips with respect to the plane of the central 
dise is always small, vanishing when w;/a is equal to 
+o, and the disturbing frequency is never 
Thus, the 


I/y or 
near the dynamic blade natural frequency. 
blade tip masses behave almost as if the blades were 
rigid. Hence, blade flexibility has little influence on 
the whirling velocities along the curves A; and Xo. 

The A; and 4 curves are quite different from their 
rigid disc counterparts for a substantial range of values 
of w;/a and, in particular, Fig. 4a shows that this dif 
ference is closely connected with blade resonance. To 
understand the influence of the blades, consider again 
the rigid disc system. The shaft exerts a moment on 
the disc during whirling which is perpendicular to the 
plane OZO’. On curves A; and \y of Fig. 3a, the mo 
ment always tends to increase x algebraically, where 
x < 0; it is always greater numerically than MM)’ = 
6B,/L,°. For the bladed disc system, the same mo 
ment exists, but because of blade motion it can go be 
low My’ numerically, although it still has the same 
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Whirling velocities of bladed disc system: « = 0, b = 
wp/a = 4. 
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sense as for the rigid disc. To the right of points B and 
C on d3 and \y, respectively, of Fig. 4a, the shaft mo- 
ment is greater than J/,)’, the blades moving in phase 
with the central disc. At points B and C blade reso- 
nance occurs, and the disc behaves as if it had infinite 
diametral inertia, the moment due to, the shaft being 
My’. The blades move 180° out of phase with the 
central disc to the left of B and C; here, the shaft mo- 
ment falls below 1)’, which means that x becomes 
positive. Thus, to the left of B and C, the blade 
motion causes the bladed disc to behave in much the 
same way as it did on the upper portion of the \» curve; 
this explains why the A; and ), curves fall below the 
comparable rigid disc curves. 

Figs. 5 and 6 give the whirling velocities as functions 
of w;/a@ for discs carrying blades having w, = 2a and 
Ww, = 4a, respectively. As the blade stiffness is in- 
creased, the departures of the \j, A», A3, Ay, curves from 
the corresponding curves for the rigid disc system de- 
crease; whatever departures do occur are easily traced 


to blade flexibility. Again, in the neighborhood of the 
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Fic. 7a. Whirling velocities of bladed disc system: ¢€ = 90°, 


O, wa/a = l. 
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bladed disc 


O, wp/a = 1. 


Mode of motion of 


b= 


Fic. 7b. system: ¢« = 9, 


points B and C, large blade amplitudes must be e) 
pected. The ratio L,x/pis found from Fig. 4b. 


Case (ii): « = 90° 
Eqs. (5) now become 


9 Q 


o102 —-6=0, Lix/p = 3/a2 


where y again equals zero and 


— |] 


| + (1/S) (1 — 2n)])d?) 
2— (J 


— 2n)d? f 


Go, = 


02 = 


In Fig. 7a, the heavy solid curves represent the values 
of X and w;/a which satisfy the first of Eqs. (9) when 
Wp = a; the light solid curves are the same as those 


given in Fig. 3a; the heavy dashed curves give the 
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values of X and w;/a at which blade resonance occurs; 
and the curves and points have been labeled as in Fig. 
4g. The mode shape, as a function of X, is given in 
Fig. 7b. Figs. 8 information as 


shown in Fig. 7a but with wy = 2a and wz = 


and 9 present 


la, re- 
spectively. 

The influence of the blades on the whirling velocities 
is sufficiently similar in the present circumstances to 
that of Case (i) so that a detailed comparison between 
it and the rigid system need not be repeated. How- 
ever, a comparison of Fig. 7a with Fig. 4a reveals a 
number of interesting differences and similarities. 

The most obvious and probably the most significant 
difference is the substantial lowering of all the \ curves, 
with the exception of dj, and the curves representing 
blade resonance conditions. The lowering of the dashed 
curves was to be expected, since this follows from the 
known result that the stiffening effect of rotation is 
more pronounced when the blade moves in a direction 
perpendicular to the plane of rotation than when the 
motion takes place in that plane. The lowering of the 
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Fic. 9. Whirling velocities of bladed disc system: « = 90°, 


b= O, wa/a = 


\ curves indicates that the effect of blade motion on 
the whirling angular velocities is appreciably greater 
in the present case than in Case (i). Another differ- 
ence is the substantial interval over which the A; and 
Ay curves are in close proximity to blade resonance con- 
ditions. This means that large blade motion, rela- 
tive to the central disc, is to be expected over a wide 
range of whirling conditions, and this, in turn, points 
to the possibility of blade fatigue under a variety of 
operating conditions. A final difference to be observed 
is that the asymptotes are approached by the dz and As; 
curves as w;—> © ; here, the asymptote has the value 2 X 
At first 


sight this appears unusual, since the dynamic natural 


V 2/73 instead of 2 as in the previous case. 
frequency of each blade rises with increasing w;, and 
this seems to imply that the bladed disc should behave 
as a rigid disc for large w;. However, the explanation 
is provided by remarking that, since p and ¢ remain 
finite as w; — ©, the magnitude of the blade disturb- 
ance remains finite while the disturbing frequency in- 
creases with w; and in the limit the ratio of the disturb 
ing frequency to the dynamic natural frequency ap- 
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Effect of change of blade angle on whirling velocity 
and mode: w; = ¢g, b = 0, wR = a. 


Fic. 10. 


proaches 1/+/2 [see Eq. (A-9)]. Thus, in the limit 
there is appreciable blade tip motion, this causing the 
departure from both Case (i) and the rigid disc system. 

The A; curve and the \4 curve when sufficiently far 
above point C are similar to their counterparts in Fig. 
4a. This is due to the fact that the effect of blade flexi- 
bility on whirling is relatively unimportant along these 
two curves. 

Fig. 7b provides a picture of the mode shape. The 
similarities and differences between the curves on this 
figure and those in Fig. 4b are readily observed and 
explained with the aid of the previous remarks. 

Figs. 8 and 9 give the whirling velocities for discs 
with blades having w, = 2a and wz = 4a, respectively. 


Case (iii): n= 1/2 


Here, since | — 2 = 0, the first of Eqs. (5) is satisfied 
for all « and \ and the remainder reduce to 


AN=1, Lix/p = 3/2, siny = 0 (11) 


The third of Eqs. (11) indicates that the shaft and disc 
axis are in the same rotating plane during whirling, and 
the first two show that the bladed disc behaves just as 
if it were a point-mass, the blade motion having no 
influence. 


Case (iv): n = 1 


The first of Eqs. (5) is satisfied for all «and \; and the 
second, third, and fourth take the forms 


o102 — 6 = (3/32)032 
lix p= 3 02 (12 
oi{sin y — (03/8) cos y| = s] 


Fig. 10a shows the values of X satisfying the first of 
Eqs. (12) asa function of e«. For all practical purposes, 
\ is independent of «. The ratio Lix/p and the angle 
y are given in Figs. 10b and 10c, respectively. Except 
for « = 0, 7/2, or, the angle y + 0, which means that 
neither the shaft nor the central disc axis lies in the 
plane OZO’. 
affected by «. 

As noted at the beginning of this section, whirling 


We also note that the ratio L;x /p is little 


can take place if ¢ equals either zero or 1/2 or if n equals 


either 0, '/2, or 1. Steady precession can take place 


therefore under no other conditions. Since most 
bladed discs are mounted on shafts with symmetrical 
elastic properties and have e different from either 0 or 
7/2, the remark just made implies that whirling cannot 
take place in such systems—1.e., the disc motion as in- 
fluenced by the flexibility of the blades is more com- 
plicated than that occurring in whirling. The analysis 
presented in this paper does not apply to this more 
general case. However, since some understanding of 
the vibrating behavior of a bladed disc with ¢« = 0), 
7/2 is of considerable technical importance, an investi- 
gation of this point is in progress. 


(4) BLApED Disc ON CANTILEVER SHAFT WITH 
5b #0 


Eqs. (A-7) and (A-S) determine A, 11x p, y, and 7 as 
functions of n, b, and the disc constants. Hence, it is 
possible to obtain, for a given disc and shaft, the whirl 
ing angular velocities and mode shapes as a function o/ 
the disc angular velocity. There are six whirling ve- 
locities for each w; and, in general, y ~ 0, r # 0. How- 
ever, since a rotating shaft with unsymmetrical elastic 
properties is rarely used in gas turbine and axial-flow 
compressor practice, this case is of little practical im- 
portance and numerical computations were not made. 


APPENDIX 


A momentary position of a moving system, col- 
sisting of a bladed disc, D, mounted on a cantilever 
shaft, S, is illustrated in Fig. 11. 

The axes OXYZ are fixed in space, and O is at the 
built-in end of the shaft. When the system is at rest, 
the disc c.g., O’, the dise axis O’P’, and the shaft axis 
lie on the line OZ. The O’X’Y’Z’ axes remain parallel 
to OX YZ during motion. We now admit the possibility 
that the planes OZO’ and O’Z’P’ may not coincide; the 
angle between these planes is designated by y. Thus, 
the shaft axis may be deformed into a space curve 
The unbent shaft is assumed cylindrical in form, and 
the plane OZQ contains one of the principal axes of 11 
ertia of the cross section when the axis of S coincides 
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with OZ, the angle between OZO’ and OZQ being r. 
The disc consists of a rigid central portion D,, around 
whose periphery is mounted a large number of evenly 
spaced blades. The central portion has mass .V and 
diametral moment of inertia /;. 

The blades are taken to consist of identical massless 
cantilever beams carrying concentrated tip masses, the 
total amount of the tip masses being 17. We also as- 
sumne that the blade cantilever consists of a sufficiently 
wide strip of thin material so that the tip mass can be 
taken to move only in the direction perpendicular to 
the plane of the unbent strip. The angle between this 
plane and the plane of D, is designated by «. Each 
blade is mounted on the central disc in such a manner 
that the unbent ax's of the blade coincides with a radial 
line through O’. Fig. 12 shows a section of the central 
disc with just one blade attached. 


We now assume that the system moves with 


2, X, ¥, T ¢, and w; = constant (A-1) 


The force F and couple N required to be exerted at O’ 
in order to have the disc move in accordance with Eqs. 
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Instantaneous position of bladed disc system. 
angle Z’O’P’ should be labeled x. 
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Fic. 12. Portion of bladed disc. 


(1) lie in the O’X’Y’ plane and are given in complex 
form by the formulas 


F —[(C, + BC,’ + i(C. — AC,’)eJe’* | 
(A-2) 
N —{—D, + BD; — i(D, + AD;)e'Je’* \ 
where 
C = (M/2) (1 + cos” €) + M,\pe? 
C, = (ML/4) (203; — ¢)¢x sin 2e 
oii / Ws), 
C;’ = Mw; \ 8 + 9 cos" e+ 
L t) 
sill € 
I{1 — (tanh o/c)]9 
D, = MLp¢? sin 2€/4 
— ML?*(2w; — ¢)¢x sin? « ‘ 
| 
lo(20s — @)%X\ (4.3) 
MLw;* cos € | f, 
D, = - 
2 Vfl — (tanh o/c) 
1 — sin? € 
| L(2w; — ¢)¢x COS «€ 
w3" — sin — 
_/[1 — (tanh a/o) : 
pe” sin € 
B = 
< L al | 
ws"  — sin-e|— B? 
J{1 — (tanh o/a)| j 





These equations were derived using the usual assump- 
tions of small motion. In general, F and N are not 
in and perpendicular to the OZO’ plane, respectively, 
and they rotate with angular velocity ¢ 
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The force and couple that the shaft exerts at O’ also Equating Eqs. (A-4) to the expressions on the right. 
lie in the O’X’Y’ plane and are given by the expres- hand side of Eqs. (A-2) and using the abbreviations 


sions 
M sin? ¢ 
12B : o=4-—-]1+ OM (1 — 2n)@ | d? 
— (1 be *'”) — aiblo 
§ ” , 
FS 207 ML? cos” € 
6B as/. -_ o= J2— od — 4 
L? (1 + be Yr") xe" | ef "9 OL)? (1 v 2p »). 1 — 2n)d 
a 
(A-4) V I 
6B 9 _= Z <3 ®) S as ») 
—4 | . - (1+ be~7#7) o- o3 ([) (;) br” sin e 
L,? 
4B wt bn a a —_— (1 — 2n)/n? A-3 
L, (1 + be ) xe e€ L : 


— (” — -) 
— sin? e — 
/[1 — (tanh a/c) ] n 


a be m 3B 


respectively, where 


B= B, + By — B, — B, og = n\ V3 er = 
- 9 ' dees s+ o Wp l MoL; 
° ° ° ie 8 A = p/ ’ yy = D 
and B, and By, are the bending stiffness of the shaft in, on ws 
and perpendicular to, respectively, the OZQ plane. m = real number 


we obtain the equations of motion 


2Qir L Pe 96(od~) . Lix 
o, + 4be~2** = 1 = o3(1 — 2n)? + 2(1 + be" +) | 6 | 


| p 
A-6) 
3 9: 2tr ° 2i(r7+y¥) 14 11x 
o3 + 31(1 + be") = t[oo + 2be "7 *] e&”? 
8 p 
In real form, these equations are equivalent to 
Lix oa, cos y + 4b cos (y + 27) 
p 2(1 + bcos 2(y + 7)] 
3 cos y + (3/8) o3(1 — 2n)? sin y + 3b cos (y + 27) 7 
_ A-( 
a2 + 2b cos 2(y + 7) 
(3/4) o3[4n(n — 1) — b(1 — 2n)? cos 27 + bcos 2(y + 7)] + 26? sin 2 (y + 27) + 2b(202 — 3) sin 27+ 
2b(o, — 3) sin2(y + 7) = 0 
o\02 — 6 — (3/32)037(1 — 22)? + (3b/4)o3[(1 — 2”)? sin 27 + sin 2(y + 7) + 2b? cos 2(y + 27)] + (As 
2b(202 — 3) cos 27 + 2b(0; — 3) cos2(y + 7) = 0 ; 
oi(sin y — bsin 27) + [(1/8)o103(1 — 2”)? — 467] cos y + (6/2)03(1 — 2n)?[cos(y + 27) + sin (y + 27 
= () 
The displacement of the blade tip mass from its equilibrium position is given by the formula 
mw = Acos(@+ yy) + Bsin(y +60+4+ y) A-9 


where 6 defines the location of the blade with respect to some radial line fixed in the central disc and y is the third 


Eulerian angle of the disc which is equal to (n — 1)¢¢. 
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Prediction of Yawing Stability Characteristics 
of Airplanes During Catapulting’ 


HENRY J. KELLEY?t 


Grumman Aircraft Engineering Corporation 


SUMMARY 


This paper presents the results of a theoretical investigation 
of the yawing behavior of an airplane during catapulting. The 
equations of motion of the airplane are developed, linearized, and 
solved approximately under reasonable physical assumptions 
The analysis is based in part on previous work by Stevens." ? 

For the relatively simple case where the tail (or nose) wheel is 
permitted to caster freely, the equations of motion are solved 
explicitly in terms of power series developments. In the more 
complicated (nonlinear) case involving a locked and skidding 
tail (or nose) wheel, the equations are solved by means of a finite 
difference procedure. Under further simplifying assumptions a 
stability criterion is derived whereby a qualitative prediction 
of the stability of the airplane can be made, and trends as param 
eters are varied can be examined. A significant effect is that, as 
the tires become more flexible, the stability of the system im 
proves. 

Although the major portion of the paper is concerned with the 
single-pendant towing arrangement, the case of the V-bridle has 
also been considered in some detail. A nonslipping bridle results 
in an extremely stable system, but an unfortunate arrangement 
combined with large misalignments can result in slipping and 
possible instability. A finite difference solution is necessary in 
the analysis of such a case. 

A comparison is shown between calculations and a dynamically 


similar model typifying current jet-fighter parameters. 


SYMBOLS 


} = forward displacement of the airplane along the 


catapult (see Fig. 1) 


X = lateral displacement of the airplane from the 
centerline of the catapult (see Fig. 1) 

V = yaw angle of airplane (see Fig. 1) 

b = bridle angle (see Fig. 1) 

= bridle length (see Fig. 1) 

h = distance from bridle connection aft to c.g. (see 
Fig. 1) 

( = distance from c.g. aft to contact point of main 
wheels (see Fig. 1) 

F = bridle tension force - 

7 = thrust force 

VW = W/g = mass of airplane 

g = acceleration of gravity 

P, = transverse force on tail wheel 

Ps, P = transverse forces on main wheels 

ly = yawing moment of inertia of airplane about c.g 


distance from main wheels aft to tail wheel 


(see Fig. 1) 
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Xo, Yo, 
WY, etc. = zero order solutions describing basic forward 
a motion of airplane on catapult 
y, etc = first-order solutions describing small disturb 
ances from zero-order solution 
a = tire slip angle due to lateral flexibility 
Q = lateral load on both tires perpendicular to 
undeformed plane of tire 
C7 = lateral flexibility coefficient of tires; radians 
slip angle per pound lateral force 
C.P = cornering power of tire (defined in reference 3 
as pounds lateral force per degree slip angle) 
d = exponent factor in assumed solution of linear 
system 
A, As, ete = coefficients of characteristic polynomial 
(Routh’s discriminant ) 
p,k = radii of gyration of airplane in yaw about point 
0 and about c.g., respectively 
t = time 
B,, Be = equation coefficients in power series method 
An, 0, = coefficients of /" in power series developments 
of y and ¢, respectively 
Yo, $0, 
Xo, vor 
etc = initial values of y, ¢, x, ¥, etc 
V = air speed due to wind and carrier forward mo 


tion 
coefficients of aerodynamic force and moment 


terms 
S = wing area 
B = wing span 
p = air density 
bus 4 = aerodynamic sideforce and yawing moment 
coefficients, positive for forces to the right 
and moments counterclockwise, respectively 
o = correction increment accounting for damping 
contributions of wing and fuselage 
A = time increment 
Ci, C2, 
D,, D:, 
Ei, Es, 
etc. = matrix elements in finite difference method 
m = coefficient of friction 
F,, F: tension forces in left and right leg of bridle, 
respectively 
B, Bi = bridle semiapex angle in general and neutral 
positions, respectively (see Fig. 2) 
R = sideforce due to V-bridle 


bridle leg angles (see Fig. 2). 
differentiation with respect to time ¢ 


¥2 = 
Dots represent 


Y 


INTRODUCTION 


y I {iE PROBLEM OF INSTABILITY IN YAW of an airplane 
during catapulting is one that has arisen in prac- 
tice only in the fairly recent past and has received com- 


) 
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Fic. 1. Single pendant arrangement. 


paratively little analytical attention in the literature 
thus far. The only previous theoretical work known to 
the writer is that of Stevens,' * of Chance Vought Air- 
craft Division. Reference 1 contains an excellent 
account of experiences with several airplanes and quali- 
tative explanations of their behavior, along with simpli- 
fied analysis (constant forward velocity, rigid tires) 
of both the single-pendant and the slipping V-bridle 
cases. In reference 2 the equations of motion are for- 
mulated to include the effects of tire flexibility and, in 
the V-bridle case, the effects of towing cable elasticity. 

It is the purpose of the present report to present 
methods for quantitative prediction of airplane behavior 
following initial misalignment on the catapult in terms 
of time histories of the motion. The equations of mo- 
tion used herein for the single-pendant case are in 
agreement with those of references 1 and 2 except with 
regard to minor tire force effects due to path curvature. 
The slipping V-bridle analysis can be shown to be in 
satisfactory correlation with that of reference 1 for 
small apex angles if certain differences in assumptions 
are accounted for. 


THEORETICAL ANALYSIS 


Equations of Motion 


Consider an airplane being towed on a catapult by 
means of a single-pendant towing bridle as shown 
schematically in plan view by Fig. 1. The point 0 is 
the intersection of a line joining the contact points of 


the main wheels with the plane of symmetry of the 
airplane. The displacement of the point 0 forwarq 
in the direction of the catapult is Y and that laterally 
to the right is X. The angle of yaw of the airplane 
measured clockwise from the direction of alignment js 
VY, and the corresponding angle made by the towing 
bridle is @. The following purely geometrical relation 
between ®, V, and_X exists: 


asin ® = X + (6+ c)sin¥ | 


Here, a is bridle length, 5 is the distance from the bridle 
connection aft to the c.g., and c is the distance from the 
c.g. aft to 0. 

Forces acting on the airplane in the Y and X diree 
tions, respectively, are given by 


—-M(Y - We sin VW — Wc cos VW) + Fos & + 
T cos ¥ — (P, + P2 + P3) sin ¥ =0 (2 


—~M(X +4- Wc cos V — Wc sin V) — F sin & + 
T sin WV + (P; + Po + P3) cos¥ =0 (3 


The mass of the airplane is 1/7. 7 is engine thrust and 
F is tension in the bridle. The forces P are reactions 
of the main wheels and tail wheel in directions per- 
pendicular to their planes. 

Moments about the point 0 are expressed in the fol- 
lowing equation: 
—(Jy + Mc?) — Mc(X cos ¥ — XW sin ¥) — 


F(b + c) sin (® + VW) + Mc( J) sin VW + 
VYwWcos¥)— Pl =0 (4 

















Fic.2. \V-bridle arrangement. 
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YAWING 


The fundamental solution of these equations describ- 
ing the desired rectilinear motion of the airplane on the 


catapult is given by 
Vo -_ 0, Po 
MY, = + T, 


Since it is desired to study the stability of the airplane 
under conditions in which this basic motion is dis- 
turbed slightly, it is permissible to assume solutions of 


= 0, Xo = 0 


Yo = Y)(t) 


the form 
Y a = Xo + x, 


® = Dy + g, 


Yo + Jy, 
H+, 


in which the variations y, x, y, etc., from the zero-order 
solutions Yo, Xo, Vo, ete., are assumed to be small of the 
first order. This procedure results in a set of linear 
differential equations if products of first-order quanti- 


ll 


Vv = F=F,+f 


ties are ignored. 


ag=x+(b+oy (5) 
— My +f — (P; + P, + P3)yp = ( (6) 
-M(¢+ cl) — he + Ty + 
(P, + P, + P3) = 0 (7) 
—(Iy + Mc*)¥ — Mcx® — Fo(b+ 0c) X 
(¢ + y) + MYocy +t. McV wo = P,l = (0 (S) 


Simplifying Assumptions 

Since the displacement y does not occur in the dynam- 
ical equations {[Eqs. (7) and (8)] governing x and y, 
it is possible to discard Eq. (6) for fore-aft motion y as 
irrelevant to the stability problem at hand. 

Further, the nature of the tire lateral forces on the 
main wheels is such as to allow the replacement of Eq. 
(7) by a constraint equation. This follows from the 
assumption that the main wheels do not skid laterally 
for motions of the magnitude under consideration. 
With rigid tires this amounts to prescribing that the 
angle of the tangent to the airplane’s path %/Y, shall 
be the same as the yaw angle y—+.e., that the airplane 
travels in the direction in which it is pointed. Although 
this assumption is made in the analysis of reference |, 
it is pointed out there that the inclusion of lateral tire 
flexibility is essential to realistic results. This will 
also be shown to be the case in the present analysis. 
The assumption then is that the path angle of the air- 
plane is equal to the yawing angle plus a slip angle due 
to tire flexibility, 


Tire lateral load and moment due to path curvature 


K=> 
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Fic. 3. Typical tire loading curve. 
%/Vo=Wta (9) 


The variation of the angle a is known to be nearly 
linear with lateral load up to the load at which skidding 
begins (see Fig. 3), according to tests described in 


reference 3. The equation 


a = ¢rQ = cr(—P2 — P3) (10) 


defines the tire flexibility coefficient cr. Here, Q is the 
reaction on the tires perpendicular to their undeformed 
plane, equal and opposite to the load (P: + P3;). In 


the present case, 


Q = 


The flexibility coefficient c; is related to the tire stiff- 
ness parameter termed Cornering Power’ as follows: 


—Fil(¢ + y) + M Vow = M(x + cy) + P, (11) 


cr = 1/[2(57.3) (C.P.)] (12) 


(d*x dY>") = x Y,? 


(which is not given by dy/d Yo as in reference 2) are omitted here as small compared to the inertia contributions 


throughout most of the run. 


t= Vilytoal—Khle +) + MY — Me + cf) + Pil} 


The equations of motion of the system with the simplifications described above are 


(13) 








Iw + Mc(é + of) + Po(b +06) (¢ + ¥) — MVYocy + Pil = 0 (14)* 


A further convenient assumption would be that the tail (or nose) wheel is not restrained from swiveling freely 
(P; = 0), which would make the system particularly amenable to mathematical analysis. In cases where the 
motion is known to be of an oscillatory character with tail wheel locked and skidding, this assumption is certainly 
a conservative one, since the action of the skidding wheel would be of a dissipative nature. The possibility exists, 
however, that a tail or nose wheel under high vertical load could resist skidding over a sufficiently large portion of 
the run so that an initial angular misalignment could result in excessively great lateral displacement. It appears, 
then, that a realistic analysis of an airplane with tail wheel locked must include the untidy case in which the tail 
wheel alternately skids and restrains the airplane. 

The analysis following is divided into four phases. In the first, the stability of the system is examined qualita. 
tively in terms of the parameters entering the equations. The second presents a method suitable for calculation 
in the free tail-wheel case which is both self-checking and of good accuracy. The case of a locked tail wheel 
is treated in the third by means of a finite difference technique. A fourth section treats the case of a slipping 
V-bridle in a similar fashion. 


Stability at Constant Forward Velocity 

While it may seem at first paradoxical to assign a constant value to the forward velocity Yo while assuming at 
the same time a finite forward acceleration VY, necessary to maintain positive bridle tension /o, a certain justifica- 
tion for this can be found by considering a dynamically similar situation in which the inertial reactions due to 
forward acceleration are replaced by forces of a gravitational nature. The motivation for this assumption is, of 
course, in the fact that equations with constant coefficients can be examined by means of Routh’s stability criteria, 
Accordingly, with P; = 0, if x is eliminated from the equations by means of the relation (5), and solutions of the 
form e“ are assumed, the following system of algebraic equations written in matrix form results: 


YoMbcrd? -f (b + c)x aa Y(1 +. Crl ) —_ VY oMacrX? —ar — VocrFo 


M([p? — c(b + c)Jd\*? + (Fob + Tc) Macd*? + Fy(b + c) 


* The coriolis acceleration term ( Yoy term) is omitted here as small. 


Requiring the determinant of the matrix of coefficients fore also be positive for stability. Substituting into 
to vanish so that the amplitudes Y, ¢ may be different this the expressions for the A’s and simplifying, the 
from zero results in a characteristic polynomial of fourth inequality that must be satisfied by the parameters of 
degree in \ set equal to zero, the system for stability is 
A,\* + Aod® + Ask? + Ask + Az = O (16) p*[p? — c? + b? + ab + (ac/Fyc,)] X 
[(6 + c)? + ab — ac (T/F))| — (p? — c?) X 
[(6 + c)? + ab — ac(T/Fy) |? — ap*td{1 + 
(T/Fo)] + [(6 + ©)/Fer]l} > 0 (19 


in which 


> 


A, = Yo M°cra(p? —= e) 


A: = Map’ 

cee gr = id Fecr)) (17) Since the other parameters are known to a fairly high 
OoOT ° . . ° . 

le Ho+ op a(n 19. oo 

As = YoFocr\ (Fo + T) + [(6 + 0)/cr]} segiinel lou shilling. This cnallicheat normally is he 


The necessary and sufficient condition for stability of quired to be greater than some certain positive value 
the system is that all of the expressions in the following for parameters representative of current aircraft. The 


sequence shall be of the same sign relation also provides a means for establishing trends 
4 of the stability of the system as other parameters are 
Ay ; 
varied. 
Ay | | 
(AsAs — AsA,)/Az wae 
en siete deal Solution for Varying Forward Velocity 
A, < [A»?A 5/ (A2A rn A,A s)] 


Since, in practice, the question of stability in the 
a criterion that is due to Routh.’ The first three of classical sense is largely an academic one, the problem 
these are always positive for airplane and catapult remains of predicting quantitatively the behavior of 
parameters in the practical range. The last, which is an airplane during catapulting so that the degree of 
generally known as Routh’s Discriminant, must there- yawing instability may be ascertained. The attend- 
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YAWING 


ant question of acceptable limits on airplane behavior 
is likewise an open one which is not likely to be settled 
until a number of borderline cases have occurred in 
practice. 

The Eqs. (13) and (14) may be handled for arbitrary 
time variation in Y) corresponding to typical catapult 
runs by means of an analog computer [but with certain 
practical difficulties if Yo(0) = 0]. An analytical solu- 
tion suitable for computing time histories manually 
can be obtained, however, if forward acceleration is 
taken to be constant at some average value and forward 
velocity is accordingly a linear function of time. The 
equations become, with the substitutions /) = Ng and 
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The coefficients B, are given in terms of the parameters 
of the system as follows: 


l : 1 a i 
B, eg ) us C B, _ + CT 

MNegcrb Mbcr 

—ada ~——@ 
B; - By; =e . 

b MNebc, 2 

—Fe bFy — cT (24) 
Bs - B, = 

Mb M(k? — bc) 
gs ac “— Fi(b + c) 

(k? — bc) M(k? — bc) 


The type of solution suggesting itself for the equations 








Y, = Net and with P; = 0, 
t= Netly + cr[—Pole + ¥) + 


in this form is development in power series, 


vy _ > a,t", . sila 2. b,t” 


MNegw — M(%® + eb)}} (20) a q (25) 
Ty) + Mc(% + of) + Folb + ©) (e + W) — Initial displacements and velocities are given by 
MNgcy = 0 (21) 


V(0) = =a, VO) = bh =a 
Xo b+ec 
¢(0) = go = + Wo = bo, 


a a 


Eliminating x by Eq. (5), these equations take the form 


ti + Bil + Boty + Bite + Bae + Be = 0 (22) ¢(0) = @ = b 


i + By + Be + Be = 0 (26) 


(23) 


Substitution of the series (25) in Eqs. (22) and (23) yields relations between the coefficients a,, 5,, since the col- 


lected coefficient of ¢” must vanish in each equation for every 7. 


Y {(n + WI) Biraygr + (Ct + Bs Ong} OAD {a + V)ndngs + Body + (nm + 1)MBabayi + Bobrajt" = 0 (27) 
=() n 1 
po ‘(mn + 2) (n + l)dyso + Bea, + (n + 2) (n + 1) Barbrgo + Bgd, }t" = 0 (28) 
={) 
If % = vo = O (the usual case), the coefficients a,, b, vanish for n odd, and the coefficients a2, b2 are obtained by 


simultaneous solution of the two equations, 


2(1 + B,)dae + 2(B; + By)by + Boao + Bsho = UV (29) 
2d2 + 2Brb. + Beato + Bgbo = 0 (30) 
For larger n(n > 4), the values of a, and b, may be computed from the following recursion formulas: 
: — B; ( B.[By + (n — 1)Bs)) 
a = Ee An-2 + 
a{Ba[(*s — 1) + B,] — By — (2 — 1)Bs} \' By(n — 1) f 
j B;[B, + (n — 1)B3}) ) : 
i, - b,-2 (31) 
{ Bz(n — 1) f 


n > 4, even; 


; r (Jp B[(n — 1) + Bil) 
, = : ce eng Ayn» + 
ni B-\(n + B,) — Bs — (n — 1)B;} { (n — 1) f 
3. — 1) + B,) 
Ip. _ Bsl(n Ut, ») (32) 
r (n — 1) j 


n > 4, even. 
Although the convergence of the series is assured mathematically for all finite /, as the elementary tests show, 
it is necessary to retain a progressively greater number of terms for practical accuracy as / increases. For param 


eter values typical of current aircraft, the number of terms required is not prohibitive. 
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Analysis with Tail Wheel Locked 


The case in which the tail wheel is assumed locked is 
complicated by the badly nonlinear misbehavior of the 
tail-wheel side load P; because of its frictional char- 
acter. For zero transverse velocity of the tail wheel, 
this load can evidently assume any value less than the 
critical for skidding in order to constrain the motion of 
the airplane. For finite transverse velocity, the load 
assumes a constant value, the sign depending on the 
direction of skidding. While the value of the kinetic 
friction load is known to be generally less than the static 
threshold value, such a refined assumption as this is 
not warranted in the present case where only rough 
estimates of either are available. 

Analysis appears to be handled best by finite differ- 
ence techniques, since transition between intervals of 
skidding and of constraint is accomplished with less 
difficulty than with series solutions. The amount of 
labor with series solution increases tremendously in 
this case, since, in addition to requiring the use of both 
odd and even powers of ¢, it may be necessary to evalu- 
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ate several different series in order to piece together 4 
single run. 

The equations of motion valid during an interval of 
skidding are Eqs. (13) and (14), with a constant valye 
of appropriate sign assigned to P;. During the inter. 
vals in which the tail wheel does not skid, the constraint 
equation, 

z—- Yy —-W=0 (33) 
specifying zero transverse velocity for the tail wheel, 
replaces the dynamical equation for moments in yaw. 

The transition from constraint to skidding is char 
acterized by an increase in the absolute magnitude of 
the tail-wheel load beyond the threshold value. The 
vanishing of the tail-wheel transverse velocity signifies 
transition from skidding to constraint. 

It is consistent to include in such an analysis the 
aerodynamic forces and moments as well, since in the 
presence of an initial 30-knot wind their effect can be of 
magnitude comparable with that of the tail wheel 
With these additional terms the equations for skidding 
motion become 


x= Yoly + cr] —Fi(¢ + vy) + MNegy — M(z + cy) + rn(V + Yo) 2p + ro V + Vy)¥ 


rn(V + Yo)e + Pil} (34) 


Iw + Mc(t + cob) + Fh(b +0 (eo + W) — MNgcd + a(V + Y%)? + 


gQ(V+t YW —ga(V + Ye + Pl =0 (35 


where the coefficients 7, g are defined in terms of aerodynamic parameters, 


1 _ Oty Epis Oty 
n= fe » 'e= pSlpA 
2 oy 2 oy 
(36 
LO, : ( 
qm = 5 pSB oy =. pSBIpA (I + @) 


cy and c;, are the conventional coefficients of the aerodynamicist referred to the point 0. 
wing span, f is air density, and /, is the distance aft from 0 of the fin aerodynamic center. 
the difference between fin-on and fin-off (fin contribution alone). 


damping contributions. 


S is the wing area,’ B is 
Prefix A indicates 
The factor o accounts for wing and fuselage 


Basic assumptions of the finite difference approximation are that velocities and accelerations over a small in- 


terval of time can be represented as the averages of the values at the end points 
With x = &, y = s,t = nA, the system of equations govern- 


orderliness obtained by the use of matrix notation. 
‘ ing skidding motion is 


: Ci C2 C3 
5 Cs Cr Cs 
x 1 0 0 
¥ 0 1 0 


where 


There is a certain advantage in 


C4 g Cs 
Co Ss Cro 
(37 
0 0 
0 y 0 
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YAWING CHARACTERISTICS 
= 72 k?) . : 
C, = : Ek we + (c? + k*)ri(V + Yo) + cai(V + v0 | 
Mk? Cr } 0 ' : 
C, = (1/Mk*) [(c? + k*)re(V + Yo) + cge(V + Yo)] 
C; = (F/Ma) [(bc/k?) — 1] 
F (bc — k?) F+T (c2 + k?) (c? + k?) . cq : 
_ ( E *) (V + Yo)? V Y5)? 
eee ee M MR°Cr Oo 
C, = (P/M) {1 + [e(l + c)/R*]} (38) 
C, = (c/MkecrVo) + (1/MkR?) (ne + qi) (V + Yo) 
C; = (—]1 Mk?) (YoC a q2) (V + Yo) 
C= — Fb/aMk? 
2 Fb(a + 6+ c) (ric + 1) (V + Yo)? c 
ie a Mk? Mk? MRcr 
Co = —P(l + 0)/MR? 
The sign of P is taken positive when the value of the tail-wheel transverse velocity ¢ = x — Yow — WV is nega- 
tive. The displacements and velocities are given in terms of average values of their derivatives by 
fats e. En Ent 
Sn+1 Sn A § n A ; atl ; 
es =5 - 5 (39) 
Xn+1 Xn Xx Xn41 
Wn4t Vn Vn Vn4t 


If the derivatives in the last column of this equation were to be calculated from values of the variables obtained by 
means of linear extrapolation, the formula would correspond to the conventional second-order Runge-Kutta ap- 
proximation. The more exact form as given, however, happens to be tractable in the present case without trial- 
and-error techniques because of the sectionally linear character of the equations. 


























Substitution for the derivatives yields the following recurrence relation for the variables at ¢ = (m + 1)A in 
terms of the values at the beginning of the preceding interval: 
rT 6G CA CA CA | ] — ¥ CA GA GA GA|T. | 7 
4° 4 2 st” rT 2 2 - 
CeA CiA c 3A CoA C.A CrA C3A CoA * 
= at Si = = Sn+l i+ Sn Cio 
9 » 9 ? 9 » » » 
7 . : . = : . 4 + (40) 
A A 
= 0 l 0 Peas 0 l 0 Xe 0 
2 2 
A A , 
e, 3 0 | Yost | 0 ; ) ; Vr 0 ; 








The elements of the matrix on the left are evaluated at f = (m + 1)A; those of the matrix on the right at ¢ = 
nd. Determination of the variables at ¢ = (nm + 1)A requires solution of a system of linear simultaneous equa- 
Since the elements depend on , solution by matrix inversion is more laborious than an iteration scheme. 


It & 


tions. 





The corresponding treatment of the case of constrained motion yields the following recurrence relation : 


to be noted that only three linear combinations of &, s, x, and y may in general be specified initially since this sys- 


tem is of only third order. 








In one of the cases of practical interest (that where V,(0) 











0), the initial values of 


fand s may not be specified independently but must satisfy the relation £) = /so. 
1 D, ° 0 =~ bE his f 14 0 0 0 Te | 
D;A D.A D;A DA A A A 
D; — D, - _ = Sn4l D; + Ds Dy + Ds D; Ds Sn 
2 2 2 2 2 2 2 
A A (41) 
on 0 l 0 Kn+t 0 l 0 Xe 
2 2 
0 a ( . ( ] 
~% 0 I ‘ Unt ) 5 ) 7 vn q 





lhe values of the coefficients D,, are given by 
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l le % 
hut Rok: hw oe. Ane oe 
Ds; = —r(V + Y>) — ra Ds = ro (V + Y,) 
crYo ( (49) 

D,; = —(F/a) (7 +6+06)//] 

= l F(6+c) ; ages : "ae 
Ds = T + ~ F Ct+ato+rag+ Ft?)  +n(V + ¥) 

CT a 


The transverse tail-wheel load is obtained from the dynamical equation in yaw as 


[P,] = [Fy Ee 3 F,] é 
5 
oe (43 
v 
where 
Mp? (b 5) 
teu Awol, be. Es 
l l a 
b+ (44 
-) 
EX, = Fp,” (a+b+c)+(F+ we 
al l 
Eliminating the accelerations. 
[Pi] = [Ei Fe E; Ey] D, —D;s 7 t ; oe 
0 0 D; Ds D; Ds é 
Fr bea dD, DD; = D, L 
—]| D; - 
0 0 O Ds 0 O Ss (45 
DD; _— Ds, DD; = dD, 
O 0 1 O 0 O 1 O x 
q 0 0 o1tLoo exile 




















The Case of a V-Bridle 


The results of the following treatment of the forces and moments due to a V-bridle are in apparent disagreement 
with the expressions given in reference | because of differences in basic assumptions. Private correspondence with 
Stevens reveals that his expressions are linearized to small apex angles, although this is not stated in reference 1. 
Agreement is obtained for small angles if account is taken of his definition for coefficient of friction, which appar- 
ently makes it a function of apex angle. The deflection effects due to cable elasticity considered in reference 2 in 
the case of a (nonslipping) V-bridle are believed to be of only secondary importance except for impractically small 
apex angles. 

Referring to Fig. 2, the semiapex angle of the bridle is 8, which has the value ) in the neutral (isosceles) posi- 
tion. The difference in tensions in the left and right cables (F,, F:) is at most equal to the tangential friction 
force on the hook at which slipping commences. 


| Fi — Fo| < u(Fi + Fe) cos B = uF (46 
. The sideforce exerted by the bridle perpendicular to the airplane centerline is 
R = F, cos y; — F2 cos y2 (47 


The angles 7; and y2 depend on the angle (¢ + y) in first-order approximation as follows: 


Tv Tv 
“ne Bo + (¢ + W) cos? Bo, yo = a i ip + yw) cos? Bg (48 


Making use of these relations, 
R = (F, — F:) sin & — (Fi + F2) cos* Be + vY) = +uFy sin By — Fy cos* &(¢ + ) (49 


The moment (clockwise) on the airplane caused by the sideforce and by the difference in longitudinal components 
of tension in the cables is 


R(b + c) + Fi sin ya tan B) — Fo sin ya tan B = R(b +c) + uFy sin Boa + Fy sin? Bale + py) = + 
uF sin B(a + 6 +c) + Fla sin? By — (6 + c) cos? Bo] (g + W) (80) 
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The plus sign is associated with cable slippage from right to left over the shuttle hook. It is to be noted that, for 
the same bridle length a and attachment point as for a single pendant, the ‘‘static’’ stability term [(¢ + y) term] 
in the case of a V-bridle has a less negative coefficient than for the single pendant. This term can even be destab- 
ilizing (positive) for geometrical arrangements that appear reasonable. In this case the motion will become 
aperiodically divergent if at any point the value of (¢ + y) becomes large enough to initiate slipping. 

For the case where the bridle does noi slip and the tires are assumed rigid, the solution for the time histories of 
motion is particularly simple. Examination of the bridle moments developed in this case will provide a basis for 
predicting at least qualitatively the slipping behavior of the bridle. The equations of motion are 


¢ + ¥ = const. = (xo/a) + [(a +°b + c)/alwo (51) 
x= Yow (52) 
With Y = Ng(t?/2), the time histories of y and x take the following form: 
ca yer Tere (53) 
x = xX + (a + b+ che {1 — eo NEte FOTO) (54) 


This highly stable type of motion corresponds to a rapid alignment of the airplane’s path with the direction of the 
catapult. The large lateral and angular accelerations developed in the first instants produce high bridle moment 
that may exceed those available from friction. The moments on the airplane are expressed in the following: 
a sin? By — (6b + c) cos? Bo] (ge + w) + 

uF, sin (a +b +c) = 0 (55) 





—M(p° — Cyd => Mcx + MNegcy a Fy 


If the bridle is not to slip, the friction (u) term must adjust itself within the + limits to cancel the other terms. 
In the event the limits are exceeded, slipping commences. For the assumed type of motion given above, the 
friction term must satisfy the following inequality to preclude slipping: 


M(p? + c2)o lalla | Negt? | 
Fy sin Bola + 6 +c) > |\— re — 1 
cocina ' rode (a+b+c) . (a+b+c) bl 
|. & + ¢) Ngt* 


Me oNge~ X#?/24 + Fla sin? By — (b+ c) cos? Bo] k gore vo (56) 


(a +b+ cc) a 0 


Inspection of this quantity indicates that a maximum will be attained initially (¢ = 0). Also, by differentiation 


with respect to time, a maximum is found to exist for a value of ¢ given by 
t= V(a+b+4+ cc) [3(p? + c*?) — cla + 6+ c)|/[Ne(p? + c?) —cla+b+4+c)] =" (57) 


In most cases of practical interest, the absolute value of the frictional moment corresponding to the latter value 
of fis smaller than that attained initially. For extremely low values of forward acceleration (Ng ~0.2g’s), frictional 
moments due to large initial lateral displacements may cause the maximum absolute value of the quantity to occur 
att = /*. However, in view of the high accelerations encountered with present airplane catapults, it is necessary 
to consider only the absolute value of the frictional moments developed initially. 

Thus, the maximum frictional moment required to preclude slipping of the V-bridle is given by 


Xo (¢ + b + c) 
x * 1 , ( Ws | 


a a 


MNeg(p? + c*)Wo 
(a+6+¢c) 


+ Fila sin? By — (6 + c) cos? Bo] | (58) 


Consideration of representative airplane data leads to the conclusion that, for the normal range of parameters, 
bridle slipping is much more apt to be initiated by angular misalignments than by lateral displacements. The 
occurrence of bridle slipping offers no cause for alarm in itself as long as the bridle arrangement provides static 
stability and as long as the friction level is high enough to damp any oscillatory instability that may be present in 
the frictionless case. 

Calculation of a run involving a slipping V-bridle is carried out in a fashion similar to that in the case of a skidding 
tail wheel, in which the motion is divided into intervals where the constraint Eqs. (51) and (52) govern, or the 
dynamical equations govern, depending, respectively, whether the bridle is holding or slipping. The transition 
from holding to slipping is characterized by an increase of the bridle moment required for constraint beyond that 
available, while the transition from slipping back to holding again is characterized by the vanishing of the quan- 
tity (¢ + ¥). Consideration of a V-bridle case in which also the tail wheel is locked introduces the possibility of 
two types of constrained motion which cannot occur simultaneously unless ¥ = y = 0. Apparently, intervals in- 
volving different types of constraint cannot, in general, merge one into the other but must be interspersed by in- 
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tervals in which both the V-bridle slips and the tail wheel skids. 
in which both the V-bridle slips and the tail wheel skids are 
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The appropriate equations for such an interya] 


& = Voly + cr[—Fy cos? Bo(y + ¥) + (Fo + T)W — M(% + ch) + Pi F uF sin B)} (59 

Iy 0+ Mc(% + ch) — Fyla sin? By) — (b+ c) cos? By] (9 + W) — (Fo + T)ewy + 
Pil + pFy sin B& (a +b+c) =0 (6 
(aerodynamic forces omitted). The tail-wheel transverse load P; takes the sign opposite to that of (* — Vw 


— hp), and the bridle friction terms take the sign of (¢ + W) where prefixed by +, or opposite to it where prefixed by 
+. The finite difference procedure for solution of these equations is similar to that presented for the tail-wheel. 


locked case. 


DISCUSSION 


The theory contained herein describes the yawing 
behavior of an airplane during catapulting within 
certain restrictions imposed by the assumptions. As 
with any linear theory, the restriction to small ampli- 
tude motion is the dominating one; the deterioration of 
accuracy with increasing amplitude is of the same order 
as for the corresponding theory of the pendulum. In 
line with this is the limitation that sideforces developed 
must be small enough to assure that skidding of the 
main wheels does not occur as per assumption. It 
appears that, for sufficiently stable oscillations of a 
magnitude acceptable in operational practice, the side- 
forces are well within the threshold skidding values. 
For unstable oscillations of large magnitude, skidding 
of the main wheels produces a relieving effect, which 
makes the theory conservative. 

Calculation of a run by means of the series method 
in which the tail wheel is assumed freely castering or 
under negligible vertical load may be carried out by an 
experienced computer in 30 to 40 man-hours. A rapid 
direct check of this method is accomplished by substi- 
tuting the series with numerical coefficients back into 
one or both of the differential equations.* In cases 
where a heavily loaded nose or tail wheel has a pre- 
ponderant effect on the motion, the somewhat more 
cumbersome finite difference method may be employed, 
which requires approximately 60 to 70 man-hours in- 
cluding checking by a separate computer. This method 
is appropriate to the case where an airplane whose 
yawing oscillation may be well damped misbehaves 
because its heavily loaded nose wheel resists skidding 
after an initial angular misalignment. 

The state of affairs concerning basic data is a rather 
unique one as contrasted with the corresponding situ- 
ations in similar fields such as flutter analysis: Every 
‘important parameter entering the equations is known 
with good accuracy with but one notable exception 
tire flexibility. It is highly desirable that a series of 
tests be performed embracing such variables as infla- 
tion, vertical load, and size for the families of tires in 


current use. Possibly some attention should be given 


*Solution by IBM punched-card machines has been accom- 
plished in 1 hour, including checks. 


to the determination of the more basic constants occur- 
ring in the theory of reference 4. 

A few comments concerning the trends of the sta- 
bility of the system as certain parameters are varied 
may be in order. While a knowledge of such trends js 
certainly helpful in laying out a catapulting arrange- 
ment, it must be recognized that no choice of airplane 
parameters (c.g. location, yawing inertia) is usually 
open and that landing gear geometry is normally goy- 
erned by conditions other than catapulting. The tow- 
ing pendant or bridle geometry is dictated in part by 
launching attitude requirements. The upshot of this 
is that little leeway is left to improve the stability of a 
troublesome airplane short of rather drastic measures. 

One major trend is that the stability of the system 
improves with increasing forward acceleration and with 
decreasing tire stiffness, the two effects being of about 
equal importance. Not only the decay or divergence 
per cycle is favorably affected by these variations; 
equally important is the decrease in the number of 
cycles occurring in a given runout distance. Thus, for 
a fixed end speed, a short catapult with high acceler- 
ation is doubly favorable. The higher vertical wheel 
loads also more nearly bottom the tires, improving the 
stability by decreasing their stiffness. It is significant, 
and to a certain extent contrary to intuition, that the 
most adverse condition encountered is the ‘“‘low g roll- 
away,’’ which is usually among the first of the catapult- 
ing trials in a carrier suitability test program. An air- 
plane that is entirely stable under normal service condi- 
tions may be alarmingly violent at sufficiently low ac- 
celerations, particularly with power on. The effect 
of thrust is appreciably destabilizing at accelerations as 
low as 2g and is preponderantly so at accelerations of 
the order '/2g encountered in rollaway shots. 

The condition of highest yawing inertia is also most 
destabilizing to the oscillatory motion. Hence, the 
maximum take-off gross weight with outboard external 
stores is the condition to be considered for analysis. 
It may be said that the shortest possible pendant at- 
taching to the airplane as far forward as can be man- 
aged provides the most stable arrangement, but this 
unfortunately conflicts with the requirements of high- 
attitude angle and stability in pitch so that this limit 


may in many cases be severe. The most effective al- 
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YAWING CHARACTERISTICS 
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though drastic fix in such cases may be a V-bridle ar- 
rangement preferably attaching forward of the c.g. 

The effect of aerodynamic forces on the damping of 
the yaw oscillation is usually small, even in the presence 
of an initial 30-knot wind. The effect of a locked and 
skidding tail wheel can be anything from negligible to 
profound, depending on the vertical load. 

A comparison of theory and experiment is shown in 
Fig. 4. The dark curve is from a calculation typifying 
current jet-fighter parameters, except possibly for a 
tire stiffness higher than typical. The curves falling 
within the cross-hatched band are test results of a dy- 
namically similar model *’ corresponding to the same 
parameters and represent various initial yaw angles 
both left and right. Plots of the parameter x/ Wo 
plotted here would fall one over the other if the system 
were truly linear. The model results were not per- 
fectly repeatable, and a certain lack of symmetry be- 
tween left and right runs persisted in spite of concerted 
efforts to align the model. The small difference be- 
tween theory and the model may be attributed to slight 
discrepancies in basic parameters (tire characteristics, 
particularly, since measurements of model tire stiffness 
were crude) and to minor sources of dissipation not 
accounted for in the theory (friction, wheel drag loads). 
The theoretical prediction is slightly conservative. 


DURING CATAPULTING 


CONCLUSIONS AND RECOMMENDATIONS 


Since the danger of instability during catapulting is 
a very real one, as evidenced by certain recent ex- 
amples,' and, since the delay brought about by an ex- 
perimental approach to the problem in case of trouble 
with an airplane prototype may seriously jeopardize 
its acceptance, it is obvious that any reasonable amount 
of preliminary calculation is well justified. The labor 
involved in qualitative prediction by means of the 
simplified criterion presented herein is a matter only of 
minutes. The calculation of a time history entails 
labor comparable to that required for a three degree of 
freedom flutter case. Setting up of the problem on an 
analog computer is probably worth while if a large 
number of time histories is required. 

It is recommended that the method of analysis pre- 
sented here be considered for general use pending full- 
scale experimental checks, since it appears that the 
assuinptions make the prediction conservative. Two 
main avenues for further investigation are open. A 
comprehensive series of tire tests yielding lateral flexi- 
bility characteristics as functions of inflation, vertical 
load, and tire size is much to be desired. A systematic 
experimental check of the theory on a full-scale air- 
plane whose tire characteristics are known would serve 
to establish the degree of reliability which can be ex- 
pected. 
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Two-Dimensional Transient Motion and 
Flutter of a Wing Having Four Degrees 
of Freedom 


HERBERT REISMANN* anp GILBERT C. BEST? 


Consolidated Vultee 


SUMMARY 


The following investigation formulates the flutter problem of a 
wing-aileron-tab combination with aerodynamic overhang in 
terms of a system of six linear ordinary differential equations of 
motion with constant coefficients. The coefficients are functions 
of the airplane parameters and may ultimately be expressed as 
functions of the forward speed of the airplane. An investigation 
of the stability of the solution of this system of differential equa 
tions as a function of the forward speed will yield the flutter speed 


as well as the flutter frequency. 
INTRODUCTION 


I REFERENCE |, THEODORSEN derives the equations 
of motion for a wing-aileron combination. In 
references 2 and 3, Theodorsen and Garrick, and Kuess- 
ner and Schwarz, respectively, give the aerodynamic 
forces acting on a wing-aileron-tab combination includ- 
ing aerodynamic balance. In both cases a particular 
case of motion was assumed, and the results were found 
to be useful for the determination of the flutter speed 
i.e., the forward speed of the airplane at which a har- 
monic motiorr of the lifting and control surfaces would 
ensue. 

In 1938, von Karman and Sears gave the transient 
aerodynamic forces acting on a wing in a state of arbi- 
trary plane motion.‘ These forces are in agreement 
with those of references 1—3 for the special case of har- 
monic motion (flutter). It is the object of this investi- 
gation to derive general equations of arbitrary motion 
for the case of a wing-aileron-tab combination including 
aerodynamic balance. Such equations will permit a 
more general, comprehensive, and simplified study of 
the flutter problem. In addition, the resulting equa- 
tions will be highly amenable to solution by an electronic 
computer. 

Because of the inherent mathematical complexity of 
the subject, the underlying development makes use of 


all derivations and results given in references 1—10. 
The chief assumptions may be summarized as follows: 
(a) The flow encountered is nonviscous, incompres- 
sible, and two-dimensional. 
(b) The displacement of the thin airfoil normal 
to its mean path is small. 
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THE WAKE EFFECT 
In reference 4, it is shown that the lift (and moment 
acting on a wing in nonuniform motion consists of three 
parts: (a) the “quasi-steady” lift, (b) the ‘‘apparent- 
mass” lift, and (c) the ‘“‘wake effect.” 
The last is expressed by von Karman and Sears as 


Ly _ 7 y(é) dé 
pv 1 Ve] 
dT o(t,) 


id | 
titty I dt; 


where y(é) corresponds to Theodorsen’s 


Plo(t — t,)] dh Iy(O) B(vt) | 
— U(xo), in | 
either case expressing the vorticity distribution in the } 
wake, and I) is the quasi-steady circulation. The 
function 1 
calculated by 
lationship between the quasi-steady circulation and 


— @® is the lift deficiency function as 
Wagner in reference 10. The re 
the vorticity distribution in the wake, as given by Eq. 


(11) of reference 4, is 


+e 

- Gm ‘p i 2 
| ¥(§) \. << ag 

An expression equivalent to Eq. (1) was obtained by 

Theodorsen for the lift induced by the vorticity in the | 

wake [Eq. (VIII), reference 1} 


P = —2zmrpvbCQ 3 
where 
. ; Xo : 
2rO0C = / U’ dx» +) J 
J} Vx? — 1 
4 Xo + l ‘ 
270 = U 1X9 
" J i — | 


A comparison of Eq. (5) with Eq. (2) results in the 


identity 
Lo == 290 0 


More ver, since 


a Xo 5 LU’ dx, J 
U dx» = 250 — (4) 
os Vx"? — 1 1 Vx? — 1 


we have, by comparison of Eq. (7) with Eq. (1) and the 


use of identity (6): 
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FLUTTER OF A WING - 


A © 


| ; = U dx = 240 — [I(t) + To(0)P(vt)] (8) 
1 Vx0? l 


where 


“! dO(t,) 
I(t) = / Wh P{v(t — ty]) dt, (9) 
0 


dt, 
The combination of Eqs. (8) and (4) results in the rela- 
tionship 
OC = Q — [1(t) + Q(0) P(2) | (10) 


Eq. (10) is fundamental in the subsequent develop- 
ment. The function QC appears in the ‘‘wake effect”’ 
portions of the expressions for aerodynamic forces and 
moments and has been evaluated for the particular case 
of harmonic motion in references | and 2. By virtue of 
Eq. (10), it is now possible to evaluate this function 
and the corresponding aerodynamic forces and mo- 
ments for arbitrary changes in the quasi-steady circu- 
lation, and this concept will be applied to the case of a 
wing-aileron-tab combination including aerodynamic 


balance. 
THE EQUATIONS OF MOTION 


We choose g, (n = 1, 2, 3, 4) as the generalized co- 
ordinates, where g; and g2 denote the dimensionless verti- 
cal deflection and rotation of a wing station, respec- 
tively, g; denotes the rotation of the aileron with re- 
spect to the wing, and gq, denotes the rotation of the tab 
with respect to the aileron. 

Introducing the value of QC as given by Eq. (10) 
into the expressions for lift and moment [| Eqs. (22)—(25) 
of reference 2], we may write the generalized forces 
Z, corresponding to the generalized displacements 
g, in the following form: 


} 


Zn = > (Minn Gi + Bangi + Kings) — 
] 


i 


B,,[I(r) + Q(0)¥(r)] 


where 
(11) 


4 


O(r) = > ag + big 


*t dQO(71) 
I(r) = / Ur P(r — 1) dn 
Jo 


dry 


The 56 parameters .V/;,, Bin, Kin, Ai, 0; (4 = 1, 2, 3, 4; 
n= 1, 2, 3, 4) are functions of the geometry of the air- 
foil, the density of the air, and the constant forward 
speed of the airplane. They may be evaluated for a 
wing-aileron-tab combination including aerodynamic 
balance by equating the coefficients of the same gen- 
eralized coordinates and their derivatives in Eqs. (22) 
25) of reference 2. 

Assuming that the wing-aileron-tab combination per- 
forms damped vibrations of small amplitude about a 
Stable equilibrium position, we may write its kinetic 
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energy, dissipation function, and potential energy, 


respectively, as 


2T = 


1 ; 
> dX mga 
i=1 j=1 
q 4 
2D =  @ > BiGiGs | (12) 
i=l j=1 


4 4 
2U= LoD kugas 
i 1 j=l 
where mij .= mji, Bij = Biv Rig = Rj; The 30 param- 
eters my, Bij, Ry (4 = 1, 2, 3, 4; 7 = 1, 2, 3, 4) char- 


acterize the inertia, damping, and stiffness properties 
of the structure, respectively. 
Substitution of Eqs. (11) and (12) into Lagrange’s 


equations 


d (oT OD ov r 
( + = Z, (13) 
dt \Og» Ogn On 


results in a system of four linear integrodifferential equa- 


tions of motion. 


{ 
»s (Min — Min)Gi + (Bin — Bin)Gi + 
1 


i 


(K;, — RindGi = B,, (I(r) + O(0) (7 (14) 


(2 = 1, 2. 3, 4). 
It has been shown by Jones? that a good approxima- 
tion to the Wagner function may be obtained by assum- 


ing that 


P(r) = ae "+ ave *” (15) 
where a, = 0.165, a2 = 0.335, y1 = 0.0455, and yo = 
0.300. Consequently, 

“r dQ 30 
I(r) = aye ' dr, + 
Jo dr 


dQ ee ee (16) 
ave : arty 
Jf 0 dr 


Q(0) (7) = QO(O) ae + QO(0) ave 


Introducing the new variables &; (2 = 1, 2),* let 
‘dQ oe ‘ 
i, = ae ' dm + O(0)a é 17) 
o dr; 


Differentiating Eqs. (17) with respect to 7, we obtain the 
differential equations 


- - dQ 
& + yi = 1S) 


dr 


and comparison of Eqs. (17) with (16) results in the 
identity 


a. a 19) 
2 : 


I + Q(0)®= 
Substitution of Eqs. (18) and (19) into Eqs. (14) 
vields the following system of differential equations of 


motion: 
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1 ) 
z= (Min Min) Qi = (Bin Bin) Gi + ! 
i=1 
(K in a Rin)Gi — By, (& =i &) = 0 

(n = 1, 2, 3, 4). (20) 

4 

Q, >. biGi + agi — be — née = 0 

i=l 

(2 = 1, 2) 


DISCUSSION OF RESULTS 


The damped, small amplitude motion of a wing- 
aileron-tab combination (including aerodynamic over- 
hang) about a position of stable equilibrium is charac- 
terized by a system of six linear homogeneous differen- 
tial equations. The S86 parameters of Eqs. (20) can be 
evaluated for particular values of the forward speed, 
air density, and airfoil geometry. 

The stability of such a system may be investigated 
in a variety of ways. For instance, to establish the 
flutter speed directly, we may assume harmonic motion 
in all coordinates. The result of such an assumption 
will be a six- by six-element determinantal (character- 
istic) equation, the elements of which, in general, will 
be complex. The determination of the flutter speed 
and frequency from such an equation is described in de- 
tail in reference 7. 

Because of the inherent complexity of the numerical 
solution of the flutter stability determinant, machine 
computation has been applied to this problem.*® 
Eqs. (20) are particularly amenable to solution by an 
electronic computer of the REAC (Reeves Electronic 
Computer) type. In the application of this device to 
the solution of Eqs. 20, each of the coefficients is ex- 
pressed by a separate potentiometer setting, making it 
possible to vary these coefficients easily. Consequently, 
if the flutter problem is formulated in such a way that 
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the coefficients of Eqs. (20) are functions of the forward 
speed of the airplane, then a series of computer soly. 
tions for a spectrum of values of the forward speed wil] 
disclose whether the ensuing motion is stable or un- 
stable for arbitrary initial conditions. In particular, 
at the flutter speed, the motion will be harmonic, and 
the frequency as obtained from the machine solution 
will be the flutter frequency. This process is explained 
in reference 8 for the special case of two degrees of free. 
dom. 
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Two Methods for Predicting Air Loads on a 
Wing in Accelerated Motion 


HOLT ASHLEY,* JOHN DUGUNDJI,?t anp DONALD O. NEILSON? 
Massachusetts Institute of Technology 


SUMMARY 


Aircraft taking off over a distance of 25 chord lengths or less 
are shown to experience an appreciable loss of lift due to unsteady 
flow effects. Two methods, the first based on Laplace-transfor- 
mation solution of the integral equation and the second a simpli- 
fied approximate iterative approach, are presented for evaluating 
the two-dimensional lift build-up during a maneuver involving 
changing forward speed and angle of attack. These schemes are 
applied to a number of examples, including Wagner’s problem, 
which permit numerical comparisons. The effect of head wind 
on the unsteady lift deficiency is discussed. 


SYMBOLS 


a = linear acceleration of airfoil along flight 
path; location of airfoil rotation axis 
a = constant reference acceleration 
A = function tabulated in reference 4 
B=R+ 
(U,2/2a9b) = dimensionless coordinate used in special case 
of constant acceleration into a headwind 
b = semichord of airfoil 
Cr = wing lift coefficient 
Cp = wing induced drag coefficient 
E(¢, k), 
F(¢, k) = incomplete elliptic integrals of the second and 


first kinds, respectively 


E(k), K(k) = complete elliptic integrals 


Ei(x) = 
x dx’ 
f ? = function tabulated in reference 7 
JOln x 
Ui 
F, = = dimensionless wind speed 
V 2ayb 
K, = modified Bessel function of the second kind 
and order n 
K = integration limit in Eq. (24) 
L=l,+4 
L, + Ly = total unsteady lift per unit span of airfoil 
L = quasi-steady-state lift 
L = lift due to noncirculatory part of the un- 
steady flow 
L:, 1 = unsteady lift correction term due to circu 
latory effects, exact and approximate 
= operator indicating the Laplace transforma- 
tion S Pe — PR: \dR 
Vy = Ty/bU* = dimensionless circulation 
p = variable in the Laplace plane; coefficient of 


t in exponential take-off law 
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R = dimensionless coordinate of downstream end 
of wake (length of take-off run in semi- 
chords) 

S = (R 4 

I t’) = distance along take-off run from starting 
point expressed in semichords 

t = time coordinate 

U = air speed 

Uo, U* = constant reference air speeds 

U; = wind speed relative to earth 

Us = constant in velocity-time relation: U = 
Ut" 

x = coordinate of points on the airfoil, positive 
aft from midchord 

A.R = wing aspect ratio 

@* = quantity measured relative to air moving 
past the airfoil when head wind is present 

a = angle of attack of airfoil 

ao = constant reference angle of attack 

¥ = “vorticity” or circulation per unit length 
chordwise on airfoil and wake 

Yo = quasi-steady-state value of ¥ 

r = total circulation about airfoil 

Io = quasi-steady-state value of T 

Ice = strength of starting vortex 

€ = asmall positive quantity 

php = 7/U* = dimensionless form of y 

g = coordinate of points in the wake, positive aft 
from midchord 

t’ = ¢/b = dimensionless wake coordinate 

p = air density, assumed constant 

® = Wagner lift-deficiency function 


INTRODUCTION 


6 bx. PROBLEM OF THE DEVELOPMENT of aerodynamic 
forces on a lifting surface during a short rapidly 
accelerated take-off run, which is recognized as im- 
portant in catapult design, took on added interest be- 
cause of the recent appearance of several aircraft in- 
tended to become air borne in distances from rest of the 
order of 100 ft. The standard method of calculating 
performance by assuming steady-state lift and induced 
drag can be definitely unconservative in its estimation 
of take-off distance in such cases. This error occurs 
because lift is reduced from its steady-state value by 
the presence of a strong wake close behind the wing, so 
that several more m.p.h. of air speed are required to 
permit the airplane safely to leave the ground. Since 
this loss of lift is often a relatively small correction, a 
simple procedure would be useful for estimating it ap- 
proximately once the general characteristics of the take- 


off run are known. In a subsequent section such a 


2 
» 
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Airfoil accelerating in two-dimensional flow; some com 
monly used symbols are illustrated 


method is presented, followed by examples of its appli- 
cation. 
completeness, however, it is felt desirable to commence 


To provide a standard of comparison and for 


with a more rigorous mathematical approach to the ac- 
celerated-airfoil problem. 


MATHEMATICAL APPROACH 


It has been shown by von Karman and Sears! that 
the transient lift on a two-dimensional flat plate airfoil 
(Fig. 1) moving in an incompressible fluid can be repre- 
sented by the sum of three lifts—namely, Lo, Li, and 
L2, where Ly is the steady-state lift, L; is the apparent 
mass contribution, and L, is the wake contribution.'; ” 
Using the notation of reference 2 but replacing the air- 
foil of semichord c by one of semichord 0, these three 


lifts can be represented as 


L(t) = pl (t)To(t) (1) 


*b 
d 
Li(t) = —p / yo(x, t) dx (2) 
dt J-» 


: (é, t) dé 
L(t) = pl oo f = (3) 
66> VP-P 


Lit) = L(t) + Lift) + L(t) (4) 


In the above equations, yo is the quasi-steady-state 
vorticity distribution over the airfoil, I) is the total 
steady-state circulation over the airfoil, and y() is the 
vorticity distribution in the wake. y(é) is related to 
the steady-state circulation through the integral equa- 


b+ Sit U at 
Tif = — ¥7(é, t) [ 
I aid 


The problem then of determining the transient lift for 


tion, 


7 686 «{5) 


Str | Ste 


+b 
r 
b 


any specified time history of forward velocity U (ft) and 
angle of attack a(t) consists principally of evaluating 
the term L.,(t), which depends on y(é, ¢) [which, in 
turn, is determined through the integral equation, Eq. 
(5)|. The remaining terms L)(¢) and L,(¢) present no 
problem, since they can be evaluated rather readily 
from y and Io, which are well known from thin airfoil 
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theory. By means of the Laplace operator methods 
as used by Sears,” one can determine the wake lift L,(; 
for any combination of forward-velocity change and 
angle-of-attack change in terms of a Duhamel integral 
with the Wagner lift-deficiency function ®(s). These 
methods are briefly reviewed here. 

Using dimensionless parameters, 


M, = T./U*%, » = y/U* 


I Ss 
R= y| U dt 
b Jo 


and &’ = £/b, Eq. (5) transforms to 


"1+ R er 
t’ + | 
M= - (¢’ _ z? 
1 J we TN 2 ; a 


Since the vortices are assumed to remain stationary 
at the position where they were shed as the airfoil 
moves ahead, the nondimensional vortex strength u is 
a function of the nondimensional distance s from the end 
Making the substitution s = 1 + R - 
t’, the above equation reduces to 


ul i OO ts Soa 3 
vig =: - q § 6 
0 mS) 'N (R — s) : 


Similarly, the wake lift Z. can be reduced to the form 


Pa 
L, = pUbU* J u(s) X 
J0 


of the wake. 


l 
V2(R — s) + (R-—s)? 


ds (7 


Both of these two equations are readily recognized as 
containing convolution integrals and, as such, are par- 
ticularly adapted to operator methods. 

Taking the Laplace transform of both equations, one 
arrives at 


2+R 
£} Mo} = —Lfyu} £ \ - 
) Ly { ta j { 
= £ > 
lpUbU*s ee on a ae 


Eliminating the » from these two equations results in 


le 9 t 1 VY 2R + Ry 
= —£ iM)! 8 
esx l 

)* (2 + R Ry 


7 
~ pUbU* 


It will be convenient to rewrite the above as 


py fe \ 
pUbu*s 


. . jd Mo\ 
(dR | 


- Mo(0) | x 
ell/V2R+ Rs « 


pLiv (2+R R\ 
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Attention is now focused on the solution of Wagner's 


problem‘ that is, the case of a sudden angle-of-attack 


change a while the velocity is a constant ly. Under 
these conditions, 
2r aol 0 dM, Oot 
M)(0) = : ‘ =(0, R= 
U* dR h 
and Eq. (9) reduces to 
) Le ) ge {1/V2R + R?} 
) 


is l2mpUs2ba0f 7 pL iy (2 aa R)/R} 


The left-hand side of the above equation can be 
recognized immediately as the Laplace transform of the 
familiar &(R The right-hand side can be 
expressed in terms of modified Bessel functions of the 


function. 
_— kind ies 
second kind to give 


igi = © (1/V2R + R®} K0(p) 
b ?P; = = rn 
pe {VW (2+ R)/R} Pl Kolp) + Kilp)} 

(11) 


This equation can be solved for #(R) by finding the 
inverse Laplace transform of the right-hand side.t 
For the purpose here, it is sufficient to substitute Eq. 
(11) into Eq. (9) to yield 


“gr, 
—plit ) — R’ ‘ (O)® 
ol [/ ar’ ek R’) dR’ + Tol o(R) | 


(12) 


L(t) = 


Eq. (12) is the familiar Duhamel integral, which 
could have been obtained directly from the step func- 
tion response. It is derived mathematically here to 
illustrate its validity despite a variable velocity L(t). 
The only effect of this variable velocity is to alter the 
semichord lengths traveled R in the familiar Duhamel 
integral. 


The evaluation of the steady-state circulation deriva- 
tive with respect to distance can be facilitated at times 
by making use of the relation 


dTo/dR = (dYo/dt)[b/U(b)| (13) 


THE APPROXIMATE METHOD 


Following a suggestion from Prof. Manfred Rauscher, 
one may simplify the mathematical manipulations of 
the preceding section by using an approximate repre- 
sentation of the wake of the airfoil. 
constant @ and variable l’, the lift Z per unit span can 


For the case of 


be written? as the sum of a circulatory term and a non- 
circulatory term, the latter depending on the instan- 
taneous acceleration a(t). 


' The analysis to here including the inversion of Eq. (11) to 


yield the 4(R) function is described by Sears in reference 2 
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. ve’ 
2rpl”ba J r = dt 
Ji Wet — 1 


Ss 
. ey 
f v(§ + D) ae’ 
l Vv ¢’? = J 


When, as in the take-off problem, the wake is of finite 
extent, the upper limit &’ = © is replaced by ¢’ = 
R + 1. The integral from 1 to R + 1 then means 
that the wake is R semichord lengths long and the air- 
foil has moved a distance 6R from its starting point. 


+ rpb’aa (14) 


If it is recognized that the theoretical two-dimensional 
steady state is Ly) = 2rpl’*ba, Eq. (14) may be written 


f vé dt’ 
L 1 oer} ba 


Fi = ; - 27 (15) 
0 4 vg + D) ae’ 2U° 
1 Vt” a ] 


The method proposed for calculating L is to approxi- 
mate y by the value it would have if the lift built up in 


as 


a steady-state fashion. Thus, 


Lo = 2xpU*ba = pUTo, Tro = 2rUba (16) 
By the law of continuity of vorticity, 
y = -—dT)/ds = —2rba(dl’'/ds 


Asa 
first approximation, y is proportional to dl’/ds, the 


s signifying the distance along the take-off run. 


constant of proportionality being unimportant because 
it cancels in Eq. (14) or (15). Second and higher ap- 
proximations to y are set equal to —dI/ds, where T = 
L,/pU, and L, is the circulatory part of the lift taken 
from the result of the preceding approximation. Cal- 
culation of L/L» in each step requires only the two inte- 
These can be carried 
out analytically for a number of special cases. More 
generally, they are done by numerical or graphical 


grations indicated in Eq. (15). 


means, in which case care must be taken when handling 
the singularity of the integrand at &’ = 1. 

Application of Eq. (15) to the take-off problem would 
involve neglecting ground effect and assuming that the 
wing of finite span behaves like an airfoil in two di- 
mensions. The second objection can be partially 
overcome by taking, for Lo, the lift appropriate to the 
three-dimensional wing, since the ratio L/L» is rather 


insensitive to aspect-ratio changes.* The induced 
drag of the finite wing is altered by the decrease in 
lift. A good approximation to the induced drag coef- 
from the steady-state 


ficient Cp, can be obtained 


formula 


Ce = C;* wA.R. 

where C, is the lift coefficient based on the circulatory 
portion L, of the unsteady lift. If the acceleration 1s 
great, there is an additional drag rpb*a*a per foot of 
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span due to virtual mass effect. Cp, can also be cor- 
rected in the usual way for the presence of the ground. 
In some cases, it may turn out that the drag variation 
during take-off is substantially different from that based 
on the steady-state assumption. Velocities and ac- 
celerations can then be completely recalculated using 
the first approximation to L/Lo, and the new take-off 
run can be reanalyzed if desired. The proposed method 
may also be applied when there is a constant or variable 
head wind, but such cases are best illustrated by the 


examples that follow. 


WAGNER'S PROBLEM 


A first simple example, which permits comparison be- 
tween the approximate results and a known exact solu- 
tion, is the aforementioned case of the airfoil jumping 
suddenly from rest to constant velocity Uy or changing 
its angle of attack suddenly. This problem was first 
solved by Wagner,‘ and reference 2 gives a table of 
L/L» accurate to three decimal places. The velocity 
only changes at the instant of starting, so the first ap- 


proximation to y for Eq. (15) is 


AUGUST, 1952 
y ~ (dU/ds) = limit (U>/e) (17 
ee oO 


The wake extends from £’ = R+ 1—etot’ =R +1, 
where R is the number of semichord lengths from the 
trailing edge of the airfoil back to the point of starting, 
The first approximation to ZL is, then, 


= limit : 
Te e~O _ +1 (£’ +t 1) . 
dé 
R+1-eV¢? — ] 
L/Lo = (R + 1)/(R + 2) (18 


Despite its simplicity, Eq. (18) has the correct initial 
value of !/z and agrees up to the first power term in 1/R 
with Sears asymptotic approximation’ to Wagner's 
function for large R. After replacing R with the vari- 
able of integration R’ and recognizing that R’ = R + 
1 — &, a second approximation to y may be found 
from Eq. (18). Substitution into Eq. (15) furnishes a 
second approximation to the lift, 


In[ (1/2)(R? + 4R + 2+ VR? + 2R VR? + 6R+8)] 


l+: ; 
Le (R + 3) VR? + 2R VR? + 6R+48 
In R+4  2(R+ 4) In[(1/2)(R? + 4R + 2+ VR? + 2R VR? + 6R+8)] 
R+3 (R +3) VR? + 2RVR?4+ 6R+8 


The first and second approximations are compared with 
the exact Wagner function in Fig. 2; some numerical 
values are given in Table 1. 

The first approximation has a maximum error of 
about + 8 per cent when R = 2 and lies within 2 per 
cent for R > 20; the second approximation has a 
maximum error of about —3 per cent when R = 1!/2 








/ 
0.6 T 
05. — 4 —__t____i—_ ma 
O lO 20 30 40 50 60 
R 
Fic. 2. Ratio of actual to steady-state lift L/Lo vs. semi- 


chord lengths R moved after a sudden jump from rest to constant 
velocity (Wagner’s problem). The first and second approxima- 
tions are compared with the exact solutions. 


and lies within | per cent for R > 6. A sudden velocity 
jump is a rather severe example on which to test the 
proposed method, so these results are by no means un- 
satisfactory. The rapid convergence of the iterative 
process is demonstrated by the good accuracy of expres- 
sion (19), and some light on the physical significance 
of this convergence can be gained by comparing Eqs. 
(18) and (19) with the exact values. By overestimat- 
ing the lift and concentrating the wake at the starting 


TABLE 1 

First Second 

Semichord Exact Approxima Approxima 
Lengths, Wagner tion, tion, 
R Function Eq. (18) Eq. (19 
0 0.500 0.500 0. 500 
1/4 0.529 0. 556 0.521 
1/, 0.556 0.600 0.547 
8/4 0.580 0.636 0.562 
l 0.601 0.667 0.581 
1'/, 0.637 0.714 0.617 
2 0.669 0.750 0.650 
4 0.758 0.833 0.742 
6 0.8138 0.875 0.806 
8 0.849 0.900 0.845 
10 0.875 0.917 0.873 
15 0.915 0.941 0.913 
20 0.937 0.955 0.935 
25 0.949 0.963 0.948 
30 0.959 0.969 0.958 
40 0.970 0.977 0.969 
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int, the first approximation yields values of L/L» 


| 7 
The second 


which grow too rapidly in the beginning. 
representation of the wake is distributed too strongly 
over its first few chord lengths, and this has the effect 
of reducing the lift below its true value for low R in the 
second approximation. A third approximation might 
be slightly high but almost indistinguishable from the 
precise result beyond one or two chord lengths. 


OTHER FLOWS STARTING FROM REST 
Let an airfoil that is originally at rest be moved for- 


ward with a velocity Ll’ = U’4t", and consider first the 


more rigorous mathematical approach. The airfoil 


)) 
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rotates about an axis located a distance ab behind the 


midchord. The angle of attack is given by the expres- 
sion a = ao(l + f(t)]. Then, 


l ” nes * 
R= | U dt = —* 
b Jo (n + 1)b 


t= [b(n + R/U" +” 


or 


As can be shown using thin-airfoil theory, 


Py = 2rbaol{l + f(t)]Ugt” + 2b2ao[(1/2 — a)](df/dt) 
dVo/dt = 2wbay} Us{1 + f(t) |nt"—! + Ugt"(df/dt) + 
b[(1/2) — a](d*f/dt?)} (20) 


Upon substitution of the above into Eq. (12), one arrives at, with ¢ = ¢(R’), 


Inp baw 


Lal n (:) (n + “f . ‘ - ( ) (n + 1) ' . R’) IR! 
os )] P( — — 
n+1\R 0 ‘s R' 
b *R df : : h l = d*f l 
©(R — R’) dR’ -— = —a , 
7h & U \2 0 dtl 


From the preceding general result, the wake lift L2 for any specific motion can be obtained. 


H F b /1 df 
-@(R — R’) dR’ — ( —a]}— (0) ®(R) (21) 
U \2 dt 


For example, for the 


case of LU’ = aot (constant acceleration) at constant angle of attack ao, the wake lift becomes, simply, 


Ll —] 


*R | 
2rpU*bay 7 2% R I V R’ 


(22) 


®(R — R’) dR’ 


In the numerical evaluation of these types of Duhamel’s integrals, care must be exercised at the lower limit R = 


(.as the integrand approaches infinity. 
ponential approximations to the Wagner function 


@(R) = 0.165e-%-4" + 0.335¢ 0.800" 


Using the above approximation, the following expression will be obtained for the Duhamel integral when a 


a, U = SF gg 


” 0.0455R 
| (R’)~”" — ) @(R — R’) dR’ = 0.165e~°-9" (0.0455) ~"/" +P 4 ( 
JI 0 


Th . - . . - z 
Through a straightforward series expansion of e*, the 


K 
e-Vn + D of ge = n+ 1 Kun 
0 n 


The above series is useful provided the last term used 
is small compared with 1.0. It is useful to employ this 
series to evaluate to a point out from the origin and 
avoid the infinite integrand difficulty. Numerical 
integration from this point out to any desired value of 
K is then perfectly feasible. 

For the case of m = 1, the function 


is tabulated.t 


t See reference 7, p. 32, for values of this integral up to 


VK = 2. 


This trouble can be averted, however, if one uses one of the known ex- 


for example, 


(23) 


) 1/(n + D ef dé + 


ctr 


*0.300R ; 

0.335e—9-3® (9.300) ~"/°" 7 1) (&) i/(m + 1 e° dt 
$ s 

0 

integral can be evaluated as 

1) - n -t 

mee (24) 
2, [(2 + 1)n + 1] 1! 

(25) 


Thus, for the case of constant acceleration at con- 


stant angle of attack, 


i — 0.7736 we *4/0.0455R ' 
= _ ; e 0.0455 e! dt _ 
2rpL™aob VR 0 
0.6117 “V0300R | 
i edt (26) 
VR J0 
Lo/2rpU*aob = 1 (27) 
L1/2rpl?ab = 1/4R (28) 
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TABLE 2 
Lo + Le Lo + Le 
2rp U*bay 2rp U*bay 
q Eqs. (26) A, Eg. 
Semichords and (27) Wagner (29) 
0 0.500 0.500 0.500 
0.5 0.534 0.549 0.554 
1 0.564 0.576 0.596 
9 0.616 0.620 0.657 
8 0.780 0.777 0.808 
10 0.807 0.831 
18 0.866 0.869 0.880 


The sum of Lo/2rpU*aob and L2/2rpU*aob represents 
the same physical quantity as Wagner’s A function, 
which was determined by Wagner for the special case of 
constant acceleration’, although Wagner’s approximate 
computation of L» leads to slightly different numerical 
values than Eq. (26). 

The constant-acceleration 
proached by some catapults; 
thrust is extremely large compared with other longi- 
tudinal forces, as with missiles, aircraft designed to be 
air borne at low speed, and certain jet-assisted take-offs. 
Hence, it may sometimes be valuable to use the approxi- 
mate method in order to have a formula that will 
shorten somewhat the labor of Eq. (26). For use in 


take-off is closely ap- 
it also occurs when the 


Eq. (16), it is observed for U’ = aot that the velocity- 
distance law is U = V 2ays and 
dU jo ao 


Y ~ Nos ~ Vor +1 — ¢’) 
Substitution into Eq. (15) yields a first approximation 
in terms of complete elliptic integrals, 


L KLVR/(R + 2)] 
—_ | + 
Lo (R + 2) E[VR (R + 2) | LR 


(29) 


In Eq. (29), R again represents the number of semi- 
chord lengths from start of the take-off run, and K( ) 
and E( ) are the complete elliptic integrals of first and 
second kinds as functions of the argument 








; 





O06 + . + + + + | 
a a Se = 
O 10 20 30 40 50 60 

R 
Fic. 3. L/Lo vs. semichord lengths R for uniform acceleration 


from rest. Exact and approximate methods compared. 
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TAKEOFF WEIGHT 


O 20 40 60 80 IOC 40 

Feet 

Fic. 4. Actual and quasi-steady-state lift forces L and 
plotted vs. distance for Helioplane take-off in calm air 


i) 


k= VR/(R + 2) 


The initial value of fhe circulatory part of L/L, is '’, 
which is the case with all take-offs from rest.° The 
noncirculatory term 1/4R is initially infinite, reflecting 
the fact that Ly) = O while L has the finite value rpb?aay. 
This term makes an almost insignificant contribution 
to L/L after the first few chord lengths. 

Table 2 compares the values of Ly + Le from Eqs 
(26) and (27) with those from Eq. (29) (after removing 
the 1/4R virtual mass term) and also with the function 
A of reference 4. 

The lift development given by Eq. (29) is plotted 
vs. Rin Fig. 3. This curve and equation were actually 
the starting point of the entire investigation. They 
represent an attempt to explain a discrepancy between 
the observed and predicted take-off distances for the 
Helioplane,® a light airplane designed by Prof. Otto 
Koppen, of the M.I.T. Aeronautical Engineering De- 
partment, to land and take off on extremely short 
fields. Using the steady-state hypothesis, the sea 
level take-off run was calculated to be about 102 ft.,7 
whereas it was measured on a calm day at about 117 ft. 
All assumptions of the present section are faily well ful- 
filled, because the plane performs the maneuver at high 
thrust with main and tail gear always in contact with 
the ground until it becomes air borne at around 30 
m.p.h. At the predicted instant of take off, the Helio 
plane has moved 20 chord lengths from rest, whence Fig. 
3 gives L/L) = 0.94. To overcome this 6 per cent loss 
of lift, the air speed must be increased by about 3 per 
cent. Employing the actual acceleration at the take- 
off instant, about 8 ft. of additional run are required; 
this increment becomes 9 ft. when the (probably too 
large) virtual mass term is removed from L/L. There- 
fore, most of the observed extra 15 ft. are attributable 
to unsteady effects. Fig. 4 shows the actual variation 


+ With the exception of unsteady flow, this analysis accounted 


for take-off effects, including ground influence, as exactly 


available methods permitted 
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PREDICTING AIR 
of L and Ly vs. distance for the Helioplane. There is 
apparently no way of overcoming this loss of lift in 
still air; it depends only on the number of chord 
lengths in the take-off distance and will always exist 
as long as substantial portions of strongly circulatory 
wake remain close behind the airfoil to interfere with 
flow over it. The effect of surface wind is discussed 


below. 
One interesting additional observation in connection 


with the Helioplane is that, although take-off distance 
was lengthened, the observed take-off time was slightly 
low. The lift loss thus accounts for an appreciable re- 
duction in induced drag below its steady-state value. 


EFFECT OF HEAD WIND 


If an airfoil has been moving along at a constant 
velocity lL’; for some time and is then accelerated or if 
it is accelerated from rest in a head wind LU), the lift 
build-up will be slightly altered from the preceding 
case. Using the rigorous mathematical approach, the 
total circulation for the system can be represented as 


b+ P ip U dt 
+ APo(t) 4 me 
I + f 7 q 
m + ¢ 
ri Yow dt = 0 (30) 


In the above, Too is the steady-state circulation that 
existed prior to ¢ = 0, and the second integral is the 
starting vortex from Ij) which has traveled a large 
Since the starting vortex is 


F + / 


di + 
—é—b 


distance downstream. 
equal to I, the above relation reduces to 


“b+ f'! U dt e+b 
AT y(t) = -f Vu \; = dé (31) 


© 
s 


This equation is identical to Eq. (5) and, as such, will 


lead to a similar expression for Zz. The only difference 
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Fic. 5. L/Lo vs. semichord lengths R for uniform-acceler- 
ation take-offs with various dimensionless head wind speeds F 
Crosses represent the points at which a typical short take-off 
airplane would become air borne. 
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Fic. 6. L/Lo vs. semichord lengths R for logarithmic take-off 


law by approximate method. 


now is that the total circulation I’) is replaced by the 
incremental circulation AT». 


Let an airfoil that is originally moving with a steady velocity U; and constant angle of attack a» be moved for- 
ward with a new velocity U’ = U,; + U,t" and new angle of attack 


aq = ao[1 + f(t)] 


Then, 


[yo = 2ebaoll + f(t)][Ur + Uat"] + 2xb a0[(1/2) — a] (df/dt) 


AT) = 2rhaofU at" + (Ui + Uat") f(t) + b[(1/2) — a](df/dt)) 


d( AT») 
dt 


, ee d os df 
= 2rbayU 4 | nt” + f(t) nt + ai + i” i + UV -—@« V2 
A ( A \é at- 


i” ieee Uit Uat" * 
R = i= 
b J othe b + (n+ 1)b 
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Substitution of the above into the general formula, Eq. (12), will give 


| 1 4 f(t)nt"™—! & + r") re (, A 

nN n — 

L,(t) aay, fi’ Us a” CO. "sae 
0 


~ ®(R — R’) dR’ — 
2rpU*bay U, + Uat" U, + Uat" (R R’) dR 


l If 
U,f(0) + b (: i a) ; (0) 


( 


> = (R) (39 
U, + Uat" 82 
In evaluating the above integral, the R — ¢ relationship must be used numerically for most values of the coef- 
ficient m. However, for the case of constant acceleration (mn = 1), the time ¢ can be expressed explicitly in terms of 
the distance R, and a simple solution exists. 
Fors = I, 
t = —(U;/ao) + V 2b/ay V (U,2/2a9b) + R 


[ it 4° (! 74 

- ) -- 

he wil -, idsamle eT tales F 

== ®(R — R’) dR’ — 
[2 0 


2p U*bao 9 luz R Ms ' 
2aob ab 
| b U;? 
\ 2a a ; , \ 2aob : b /1 df 
®(R — R’) dR’ — f(0) += ( —a)— (0)} ®(R) (33 
lur . 0 dt U;? “ U, \2 dt 
re 
Voa0b \ 2aob * 
For the special case of U = U; + aot and no change in angle of attack, Eq. (33) reduces to 
le —1 r. —" 
= 7 7—, 0(B — B’) dB (34) 
2p [ 2hay 2 Vv B B-R Vv B 


where 
R = (U,t/b) + (aof?/2b), B= R+ (U,?/2aob) 


This is almost similar to the lift build-up starting from rest {Eq. (22)], except for the different limits of integration 
and the replacing of R in Eq. (22) by B. 

The constant-acceleration take-off into a head wind is also readily treated by the approximate method. In the 
first approximation of Eq. (15), circulation per unit length y is clearly proportional to dU’,/ds,,, where subscript w 
means that the quantity is measured relative to the air moving past. Both integrals must be evaluated from 
trailing edge ’ = 1 tot’ = R, +1. For take-off at constant acceleration dp, it is easily shown that 

: FE dU, dy 
U, = VU? + 2agsy, y~ = 5 
ds, VU + aS, 


This may also be written as 
/ 79 
i eas i; FY? — (Sy b) 
F, = U,/¥V 2aob being sort of Froude Number or dimensionless wind speed based on dp and chord 2b. With the 
remote starting vortex included in each integration, Eq. (15) gives the following first approximation to L /Lo: 
by F(sin do, k) l 


=1-— ; + a= 
Lo F\V F?2 + Ry + 2+ [Fit + Reo + 2)[E(k) — E(sin gy, k)] 9 4(F? + Row) 





where 


k 


V (Fi + Ry)/(Fi2 + Rw + 2) 
sin ¢, = F,/V F,;? + R, 
sin do = V(F2 + Ry + 2)R,/(Fi2 + Ry)(Re + 2) 


In Eq. (35), Ry is the number of semichord lengths of wind which has passed the wing, F(sin ¢, ) and E(sin @, k) 
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incomplete elliptic integrals of the first and second kinds, and /(k) is the complete elliptic integral of the second 


are 
kind. 


For given F;, the number of semichord lengths R of take-off run over the ground is related to R,, by 
R = Pa + 2F;? —_ 2FiV Re + F;? (36) 


In Fig. 5, L. Lo is plotted vs. R for FP; = 0, 1, 3, and 6. If an airplane of 7-ft. chord accelerates at 7 ft. per sec’ 
per sec., the corresponding wind speeds are LU’; = 0, 7, 21, and 42 ft. per sec. It is evident from Fig. 5 that a head 
wind causes the lift to build up much more rapidly over the first few chord lengths of the take-off run. However, 
because the wind also shortens take-off distance, it is not necessarily true that L/L» at the instant of becoming 
air borne is much larger than in still air. This point is illustrated by the crosses marked ‘‘airborne”’ in Fig. 5, 
which show roughly the positions where an airplane with a steady-state take-off speed of 30 m.p.h., 26 = 7 ft. and 
dy = 7 ft. per sec. per sec., would leave the ground in different head winds. Clearly, the unsteady effects are just 


about as important with a wind as without it. 


ADDITIONAL EXAMPLES 


Employing the approximate method, there are several other take-off laws that lead to simple analytic represen- 


tations of the first approximation to L//». A rather artificial exponential take-off from small initial velocity Up, 
U' = Uve”’, produces a constant value of dl’/ds. Hence, y =constant and Eq. (15) gives 

L V R? + 2R l _ 

= + (337 ) 


Ion VR?4+2R+cosh(R+1) 2R +4 (2U./pb) 


The circulatory portion of this result may be compared with the following more exact one, based on Eqs. (12) 


and (23): 


Lo + Le 4.743 (pb/ Uo) l - pb salt r _ pb ee 
—— =} — — : (0.165 3.626 “)e 0.0455 4 (0.339 1.117 “—) @-0.s000 
2xpU* aod 1 + (pb/Uo)R 1 + (pb/U»)R Uo Us 


(38) 
A related series of take-off laws in constant surface winds LU’; yields approximate expressions similar to Eq. (37). 
Thus, if the velocity of the airplane relative to the air Ll’, is given by 
Uy = Uye™ (39) 
a constant value of dU’,,/ds, is obtained. The first approximation to L/L is 
Bb , (Ui/pib) + V R,? + 2Ry l (40) 


+ : 
Io  (U,/pib) + cosh! (Ry + 1) + VR? + 2R, 2Rw + (2Ui/prid) 


where, as above, R,, is the number of semichord lengths of air having passed the wing since start of take-off. 
Aside from R,,, L/ Lo is seen to depend only on the single ‘“‘wind-acceleration parameter’ bp,/U,. The relation be- 
tween R,, and semichord lengths R traveled over the ground also involves only bp,;/ U4. 


R = R, — (U1/bp,) In [1 + (0f1/ U1) Re] (41) 
Another more realistic take-off law is defined by the following relation between U and distance s: 
U = UpIn [1 + (s/d)] (42) 


By assuming the run to start at a small initial velocity, Eq. (42) can be shown to lead to a velocity-time relation 
of the form 
Uot ii U U mi Uot (43 
= fi— or -~ = i 43) 
b Uo l 0 b 


Here, the Ai function is defined by Jahnke and Emde’ as 


= ° &’ 
Ki (In x) = > 
0 Inx 


The acceleration decreases with time in a manner that can be made to approximate closely actual take-offs by ap- 
propriate choice of the constant U. Details of the distance-time and acceleration-time relations are interesting 


mathematically but are omitted here as unessential. 
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The take-off law [Eq. (42)] leads to a first approximate L/L», 


L l 
Lo (2.90 + 2R) In (1.45 + R) si 


R 2) 1 |e +1) 
ee a ) + R2+3R+1 


R ayn [OR 1) 
ee R?+3R +1 


Fig. 6 presents a plot of L/L» vs. R for this case. 


(R + 1)VR?2 + 2RV R? + 4R + 3 


(R + 1)VR? + 2RVR? + 4R +3 


a VR?+4R+4+31n [((R + 1) + VR? +p 


i- VR? + 4R 4+ 31n [(R + 1) + V R24 9p) 


The comparatively slow acceleration is reflected by values of 


L/L» higher than in most other examples; with the virtual mass term included, the minimum is about 0.855, oc. 


curring only half a chord length after the start. 


DISCUSSION 


The exact method described in the first section can be 


applied to any acceleration and any angle-of-attack 
change of an airfoil to find the transient lift. All that 
is required is the steady-state vorticity distribution 
yo(x) for the motion considered. From this, the total 
steady-state circulation I) can be obtained, the virtual 
mass terms can also be found, and the wake lift can be 
determined according to the procedure outlined before. 
The transient lift due to various combinations of mo- 
tion can also be determined by this method. 

The calculation of the pitching moment on the airfoil 
follows simply, since it has been shown that the wake 
lift L. always acts at the quarter-chord point. The 
steady-state and virtual mass contributions to the 
moment follow readily from thin airfoil theory. 

Certain cases will arise, moreover, where it is desir- 
able to obtain a quick estimate of lift build-up during 
acceleration without resorting to the often laborious 
Approximate formulas, Eqs. (15) and 
The» are 


precise results. 
(16), have been devised for this purpose. 
by no means limited to the take-off problem or to con- 
stant angle of attack. Supposing the airfoil to be per- 
forming small changes in a, executing small linear mo- 
tions perpendicular to the mean flight direction, etc., 
one may readily calculate the pseudo-steady-state 
variation of lift and circulation Ty. Knowledge of the 
velocity régime then permits —dIy/ds to be found as a 
first approximation to y. From the result of Eq. (15), 
second and higher approximations can be obtained if 
desired. 

The examples presented in this paper show that the 
first approximation to L/L») becomes increasingly ac- 
‘curate as R— o but has a tendency to be too large by 
several per cent at lower values. This error is toler- 
able in a study of lift loss during take-off, but, if one 
desires to compute aerodynamic loading for structural 
purposes, the entire history of the forces must be known 


to the same order of precision. Such problems call for 


at least a second iteration or for application of the exact 
method. The results of first approximate calcula- 
tions will always be on the conservative side, however, 
so they are suitable for preliminary structural design in 
such situations as the lift development during a sharp 
pull-out. 

One other problem to which both proposed methods 
are well suited is the growth of aerodynamic forces fol- 
lowing passage of an airplane into a so-called “horizon- 
tal gust.” Rapid air-speed changes of large magnitude 
will occur immediately after envelopment by an intense 
blast wave produced by some sort of large explosion. 
Depending on the flight direction, air-speed increases, 
decreases, or both could take place. Accurate estimates 
of force and moment variation are necessary for control, 
stability, and structural reasons. These can be calcu 
lated preferably by at least two iterations in Eq. (15). If 
one wishes to find the angle-of-attack changes, wing 
twist, etc., after gust impact, they can be approximated 
from the result of the first iteration or determined more 
precisely by including the influence of body freedoms 
and structural flexibility in each successive approxima 
tion. 
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Simplified Laminar Boundary-Laver 
Calculations for Bodies of Revolution and 


for Yawed Wings 


NICHOLAS ROTT? anno L. F. CRABTREE? 
Cornell University 


SUMMARY 


Since the introduction of momentum methods in boundary- 
layer calculations by von Karman and Pohlhausen, many im- 
provements have been proposed. An especially simple solution 
reduces the problem to a quadrature. Here, it is proposed to ex- 
tend these methods to elementary three-dimensional cases and to 
compressible laminar boundary-layer calculations. For com- 
parison, the corresponding problems for the turbulent boundary 


layer are also discussed briefly. 


List OF SYMBOLS 


a constant, Eqs. (1.1) and (1.13) 


¢ = 

A = an arbitrary reference length, Eq. (2.7) 

b = aconstant, Eqs. (1.1) and (1.13) 

¢ = aconstant, Eq. (1.17) 

C = aconstant, Eq. (5.4) 

f = velocity profile family of the chordwise boundary layer, 
Eq. (1.5) 

F = a function of m, Eq. (3.18) 

g = velocity profile of the spanwise boundary layer, Eq. 
(3.5) 

H = ratio of displacement and momentum thicknesses of the 
boundary layer 

K = aconstant, Eq. (1.17) 

! = nondimensional shearing-stress parameter, Eqs. (1.5a) 
and (1.16) 

L = “growth parameter’ of the momentum thickness, 
Eq. (1.10) 

m = modified Pohlhausen parameter, Eq. (1.6) 

M = Mach Number 

q = 6,/0, = ratio of momentum thicknesses of the boundary 
layer in the spanwise and chordwise directions 

r = distance of a point on the surface from the axis of a body 
of revolution 

R = “growth parameter” of the momentum thickness (bodies 
of revolution), Eq. (2.3) 

S = 6,,/0; = ratio of ‘mixed’? momentum thickness to the 
“chordwise’’ momentum thickness 

T = absolute temperature 

u = chordwise velocity component inside the boundary layer 

U = chordwise velocity component outside the boundary 
layer 

? = spanwise velocity component inside the boundary layer 

V = spanwise velocity component outside the boundary layer 

* = coordinate measured along the surface in the chordwise 


direction 
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y = coordinate measured in the spanwise direction 

s = coordinate measured normal to the surface 

8 = semiwedge angle 

y = cone half-angle [in Section (2)] 

y = ratio of specific heats [in Section (4)] 

5* = displacement thickness of the boundary layer, Eq. 
(4.13) 

ey u 

0, = 1 — dz = momentum thickness of the 
f o(t- o) 
boundary layer in the chordwise direction 

é, = momentum thickness of the boundary layer in the span- 
wise direction, Eq. (3.6) 

6,, = “mixed”? momentum thickness of the boundary layer, 
Eq. (3.2) 

\ = the Pohlhausen parameter 

uw = coefficient of viscosity 

vy = kinematic viscosity 

p = density 

7 = shearing stress at the wall 


(1) THe Two-DIMENSIONAL CASI 


Sve THE INTRODUCTION of momentum methods in 
boundary-layer theory by von Karman! and Pohl- 
hausen’ in 1921, this concept has found many applica- 
tions for approximate solutions of the problem, and 
many improvements and simplifications have also been 
proposed. In the case of the two-dimensional incom- 
pressible laminar layer originally treated by von Kar- 
man and Pohlhausen, an especially simple solution was 
proposed by several authors, in which the momentum 
thickness 6, is obtained by a simple quadrature 
namely, by formulas of the type 

x 

U°—! dx (1.1) 


6,2? = avU~ 
0 
where a and / are constants, v is the kinematic viscosity, 
and Lis the prescribed velocity outside of the boundary 
layer along the body. The integration begins at the 
0, and is performed along the are 
Formulas of the type of Eq. 


stagnation point, x = 
length x on the surface. 
(1.1) have been established by Young and Winter- 
bottom’ with a = 0.47, b = 6.28; by Walz! witha = 
0.47, b = 6; and by Thwaites® with a = 0.45, b = 6.** 
We will follow here Thwaites’ discussion, who not 

** The authors’ attention was called to early Japanese publica- 
tions containing Eq. (1.1) by Hudimoto*® and Tani.*! In the 
latter paper, the values a = 0.441 and 6 = 6 were proposed, 
based on Howarth’s solution. 








554 JOURNAL OF 


only simplified the Karman-Pohlhausen method but 
discussed the possibility of certain improvements. 

The Karman-Pohlhausen method starts from the 
momentum equation, 


6, dU 
p24 4 oy, t¥ = y (*) (1.2) 
dx dx 02/.=0 


together with the boundary condition at the wall, 
v(0°u/02"),-9 = —U(dU/dx) (1.3) 


Eqs. (1.2) and (1.3) may be combined to give 


U? db, _ (s*) + 0,(H + 2) (<) (1.4) 
vy dx Oz z=0 0 z=0 


i.e., a differential equation for 6, in which three quanti- 
ties enter, which are derived from the velocity profile 
u(z)—namely, its first and second derivative at the 


wall and the ratio // between displacement and 


momentum thicknesses. The basic simplifying as- 
sumption of the Karman-Pohlhausen method is that 
any velocity profile in the laminar boundary layer may 
be represented as a member of a one-parameter family 


of profiles if it is reduced to the proper scale, or 


u/U = f[(z/6,), a] (1.5) 


where L’ and @, represent the scale of velocity and 
length, respectively, and a is a parameter not yet speci- 
fied. Suppose / in Eq. (1.5) to be expanded in the 


form 


u 


z l s\* 
- = l(a) + m(a) ( ) : ae (1.5a) 
l 6, 2! 


6, 
Introduced in Eq. (1.3), this gives 


m = —(dU/dx)(0,?/v) (1.6) 


a quantity recognized as a (modified) Pohlhausen param- 
eter. Thus, the quantity m may be chosen to charac- 
terize the members of the profile family [Eq. (1.5)] 
and will be identified with a. Because of the assump- 
tion of a one-parameter family, / and 77 are functions 


of m alone, 
l 
H = H(m) 


II 


l(m) (1.7) 
(1.8) 
functions, which may be given as soon as the family 
(1.5) is specified. 
By use of Eq. (1.5a), Eq. (4) becomes 
U , de, 
6, ‘ =]/+ (H + 2)m 


v ax 


(1.9) 
The nondimensional quantity (a ‘‘growth parameter’’ 
of the momentum thickness) 


L = (U/v)(d0,?/dx) (1.10) 


appears on the left-hand side of Eq. (9), which may be 
written in the form 
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L = 21+ 2(H + 2)m (1.1) 
or 
L = Lim) (1.19 


A universal relationship [Eq. (1.12) ] between the quan. 
tities m (Eq. (6)] and LZ [Eq. (1.10)] is a main conse 
quence of our assumptions. 

Now, using the Pohlhausen quartic’ for the profile 
family [Eq. (1.5)], the functions (1.7) and (1.8) can be 
evaluated, and by Eq. (1.11) L(m) is immediately given, 
Any one-parameter family gives this type of relation. 
ship. Moreover, it is possible to evaluate the actual 
functions [ Eqs. (1.7), (1.8),and (1.11)] for the fairly large 
number of known exact solutions of the laminar bound. 
ary-layer problem, thus checking the basic assumptions 
of the approximate method rather than the final results 
Thwaites® has carried out this discussion, and his re 
sults may be summarized as follows: 

(a) In the domain of pressure drop [negative m, Eq, 
(1.6)], the different exact solutions yield a fairly uni- 
versal relationship /(m) and //(m) and thus also L(m), 
and these functions are also well represented by the 
relations derived from the Pohlhausen-quartic. 

(b) In the pressure rise domain, the functions /(m 
and H1(m) show a much larger scatter for different cases, 
so that the dependences may be considered as universal 
only to a rough approximation. Fortunately, the de- 
pendence L(m) turns out to be of much higher uni- 
formity in the different cases than the particular fune- 
tions (1.7) and (1.8) 

(c) The relationship L(m), which has been recognized 
as fairly universal for positive and negative m, may be 
represented with good accuracy as a linear function, 

L=a-+bm (1.13 

The simple formulas of the type of Eq. (1.1) are based 
on this linearity, because introducing 1 and m in Eq 
(1.13), one obtains the equation 

U(6,2)’ + 66,2U’ = av (1.18a 
or 

(d/dx)(0,2U") = avU* 
which yields Eq. (1.1). Adopting the ‘‘best values 

for a and } from Thwaites’ discussion, 
L = 0.45 + 6m (1.14 

% 
6,2 = 0.45vU- / UL dx (1.15 
ee ( 


(d) The function /(m), whose universality for positive 
m appears questionable, determines the shearing stress 
at the wall, 


7 = ul(Ou Oz), qo = u(U 6,)1 (1.16 


and the vanishing of / fixes the laminar separation 
point. As soon as 6, is determined from Eq. (1.15), 


the value of m is known along the body, and the “uni- 
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versal” function (1.7) gives /. If the evaluation of 
i(m) is based on the Pohlhausen quartic, laminar separa- 
tion—i.e., / = O—occurs for m = 0.157. It is known 
that, for pressure distributions with a sharp pressure 
drop at the beginning, followed by a comparatively 
smooth pressure rise, the Pohlhausen method overesti- 
mates the length up to the laminar separation point. 
Thwaites, therefore, has proposed the adoption of a 
modified curve /(m), mainly based on the solution for 
Schubauer’s ellipse® 7 which gives / = 0 for m = 0.082. 
a table of /(m) may be found in Thwaites’ 
I(m) according to Pohlhausen is tabulated, 
By the use of 


[Fig. 1; 
paper, 
e.g., in Schlichting’s book,* p. 198.] 
Thwaites’ function, better results may be expected for 
pressure distributions occurring on airfoil-type bodies. 
For Hiemenz’ observed pressure distribution on a cir- 
cular cylinder, which was used as a test case by 
earlier investigators,” ' Thwaites’ function would 
give a too-early separation. This appears, however, 
to be less significant from the practical point of view, 
since for this type of flow the distance between minimum 
pressure point and laminar separation point is small 
anyhow, as compared with the total length of laminar- 
layer region. The adoption of Thwaites’ proposal 
(equivalent to the earlier suggestions to take A = —5 
or —6 instead of \ = — 12 for laminar separation in the 
Pohlhausen method) assures practical improvement. 
(Recently, the influence of the second derivative of U 
on the separation value of m has been investigated by 
Tani.*”) 

An important group of exact solutions for the discus- 
sion of the relationships /(m), H(m), and L(m) are 
those of Falkner and Skan,!! calculated and tabulated 
by Hartree,'? for velocity distributions 

U = cx® (1.17) 
which occur in incompressible flow over an infinite 


[The relationship between A and the semi- 
included special 


wedge. 
wedge angle Bis K = B/(r — 8); 
cases are K = 0 or the flat plate and A = 1 or stagna- 
tion point.| These layers are characterized by the 
fact that the velocity profiles are similar for any x, so 
that, for a certain value of A, m is a constant, as are 
|, H, and L. [From m = constant, it follows that 


6, ~ x7-*%)2) Generally, there follows from Eq. 
(1.6), 

dm/dx = (U"/U')m — (U’/U)L 
For U = cx“, dm/dx = 0 and 


(K —1)m— KL=0 (1.18) 


With the approximation L = a + bm, as in Eq. (1.13), 

this yields 

m= — nal , L= a a coe (1.19) 
(6—1)K+1 (6—1)K +1 

The same relations are obtained, naturally, by intro- 

duction of Eq. (1.17) in Eq. (1.1); this gives 
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Fic. 1. The dimensionless shearing stress / as function of the 
shape parameter m, as evaluated from the Pohlhausen quartic 
and the modification proposed by Thwaites 


U0,2/vx = a/[(b — 1)K + 1] (1.19a) 


which is equivalent to Eq. (1.19). The values in Eq. 
(1.19) may be checked against the exact results derived 
from the Falkner-Skan profiles. Actually, the check 
is a part of the basic discussion and the agreement is 
satisfactory.® 

The value of m at the stagnation point (A = 1) 
follows from Eq. (1.19): m = —a/b = —0.075 using 
Thwaites’ values for a and bd. 
beginning of the calculation; 
that, by differentiation of Eq. (13a), 


U(0,7)" + (6 + 1)U'(0,”)’ + 6U"6,? = 0 = (1.13b) 


This is important for the 
it is also useful to note 


The initial value of (@,”)’ is obtained for the stagnation 
point U = 0, as the first term drops out. If U” = 0, 
(6,2)' = 0 and a further differentiation gives (0,”)”. 

In summary, it is seen that simple and useful meth- 
ods exist for the calculation of two-dimensional in- 
compressible boundary-layer flow. In the following 
sections, analogous methods will be derived for certain 
elementary three-dimensional cases and for special 
cases of compressible flow. 


(2) Bopirs OF REVOLUTION 


For bodies of revolution, the momentum equation 
[Eq. (1.2)] must be changed to 


oe) + or + 2)0,U =» (*) (2.1) 


r dx dx 02/.=0 


where,” 1- :e distance of a point on the surface from 
the axis and x is again the arc length along the surface. 
(See, eg Schlichting,’ p. 204.) The boundary con. 
dition [Eq. (i.3)] remains unchanged. Making the 
same simplifying assumption about the existence of a 
one-parameter family of velocity profiles, the argu- 
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ments of Section (1) remain unchanged, and they lead, 


instead of Eq. (1.9), to the equation 


U 6, d(0,r) 


yr ax 


=1+ (H+ 2)m (2.2) 


This suggests the use of the following “‘growth param- 
eter’ R, instead of L: 


R = (U/vr?)[d(0,r)?/dx] (2:3) 
with this, Eq. (2.2) becomes 
R = 2/ + 2(H + 2)m (2.4) 


Now, if the one-parameter family may be taken to be 
the same as in the two-dimensional case, the function 
R(m) is equal to L(m), or, with the same approxima- 
tion as Eq. (1.13), 


R(m) =a + bm (2.5) 
i.e., with Eq. (2.3) and (1.6), 
U(6,2r?)’ + b6U'6,2r? = avr? 


or 


d (U°6,277) = avr?U"-! 
dx 


which gives, upon integration, 


6,” = avr~?U~ | r?>U"—! dx (2.6) 
0 


Thus, the equally simple equivalent to Eq. (1.1) is es- 
tablished for bodies of revolution. The remaining 
question is whether it is legitimate to use the same one- 
parameter family—t.e., the same function for R(m) 
as for L(m) in the two cases. 

Mangler!’ has given a transformation that reduces a 
boundary-layer problem for a body of revolution to 
a two-dimensional problem and which is exact within 
the limits of boundary-layer theory. Mangler’s trans- 


formation is 
, l aio _ r - 
t= — redx, 2=—2 (2.7) 
A?2 Jo A 


where A is an (arbitrary) reference length and % and 2 
are the coordinates in the two-dimensional case. Ap- 
plying Eq. (1.1) to the transformed plane flow, 


6.2 = avU~ f U’—! dz (2.8) 
0 


Here, 6, is the momentum thickness in the two-dimen- 


sional case, which is related to @, in the original axisym- 
metric flow, according to the second of Eqs. (2.7), by 

6, = (r/A)0, (2.9) 
The first of Eqs. (2.7) gives, for the differential 


d& = (r°/A*) dx (2.10) 


Introduction of Eqs. (2.9) and (2.10) in Eq. (2.8) gives 
Eq. (2.6). For the plane flow, which was obtained 
after transformation, the use of Eq. (1.1) (and, im. 
plicitly, the same family of profiles) was certainly justi. 
fied; thus, the use of Eq. (2.6) is equivalent to the 
Mangler transformation, without the necessity for 
carrying it out. 

It also follows from the Mangler transformation, but 
may also be seen directly, that the velocity distriby. 
tions [Eq. (1.17)] again give “‘similar”’ profiles along the 
surface, actually those of the Falkner and Skan family. 
The velocity distribution [Eq. (1.17) 
by an incompressible flow about a cone, but the con- 
nection between the half-angie y and the exponent K 
is different from the two-dimensional relationship. It 





iS now realized 


is, in fact, transcendental, since it involves the zeros 
of spherical harmonics of nonintegral order; it was 
calculated by Mangler and Leuteritz.'4 For small 
y, the relationship is approximately K = y?/2 (y¥ in 
radians); actually, this approximation is sufficien 
for all practical purposes, since it is known that, fo 
y = 7/2 or the flat plate, a stagnation point occurs wit 
K = 1. It is easy to draw a curve approximating the 
function AK = y?/2 for small y and going through the 
point K = 1, y = 7/2. 

A further difference between the two-dimensional 
and axisymmetric cases lies in the relationship between 
mand K. If, for a cone, r = x sin y, according to the 
first Mangler transformation [Eq. (26)], 


& = (sin? y/3A?)x' 


~ y r sak/S -: 
and Eq. (17) becomes U = ét”’”, where € is a new con- 


stant. In the transformed case, 


ww: (dU /dz)(6,27/v) = m(K/3) 


But the left-hand expression is, according to Eqs. (2.8 
and (2.9), equal to m, so that 


m(K) = m(K/3) (2.11 
This rule, applied to Eq. (1.19), gives 
m= —akK [((6 — 1)K + 3] (2.13 


which may be derived also by use of Eq. (2.6), as must 
be the case. . 

At the stagnation point, K = 1, the starting value of 
misnowm = —a/(b + 2) = —0.05625. The starting 
values for pointed bodies of revolution follow from Eq. 
(2.12) and the connection between K and y. 


(3) THe YAWED INFINITE CYLINDER 


The investigations of Prandtl," Jones,'* and Sears” 
have revealed that the laminar boundary layer on a 
yawed infinite cylinder may be determined in a simple 
way. Outside the boundary layer, in the inviscid 
region, the flow is obtained by superposition of a two- 
dimensional flow perpendicular to the cylindrical axis 


(“‘chordwise’’ flow) and a “spanwise’’ flow with a 
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LAMINAR 


velocity V constant everywhere. Again let x be the 


arc length on the surface, measured in the chordwise 
direction ; 
boundary layer is given as a function of x from the two- 
dimensional case. Let y be the coordinate in the span- 
i.e., in the direction of the cylindrical 


the chordwise component UL’ outside the 


wise direction 
axis. 

The boundary-layer problem becomes soluble after 
recognizing that any quantity, inside as well as outside 
the boundary layer, is independent of y. This argu- 
Let us con- 
then, in 


ment requires a moment's consideration. 
sider, for example, vanishing chordwise flow; 
a purely spanwise flow, we should expect a growth of 
the boundary layer in the y-direction and thus no in- 
dependence of y. Actually, we must exclude the spe- 
cial case of vanishing chordwise flow; 
flow exists, it will “blow away”’ fluid contained in the 
spanwise layer so that a constant thickness in the y- 


if the chordwise 


direction (and also an independence of any quantity 
of y) becomes possible. 

The spanwise component of the velocity in the bound- 
ary layer, v, then becomes a function of x alone. Once 
this is recognized, it may be easily concluded (see 
references 15, 16, and 17) that the chordwise flow may 
be calculated as if no spanwise flow were present. Sub- 
sequently, the calculation of v can be performed. 

For the approximate calculation of the spanwise 
boundary layer, Wild'* has proposed a ‘‘Pohlhausen- 
method,” starting from the momentum equation of 


39 


the spanwise layer,'® 


d u(V — v) dz = pv (*) (3.1) 
Oz 


dx Pe Al) z=0 


Oi. = | os (1 — 7) de (3.2) 
. V 


and the momentum equation becomes 


d (U@,,) = d (*) (3.3) 
dx ’ Oz 


z=0 


Introduce 


The boundary condition corresponding to Eq. (1.3) is, 
(O°v 08°). 36 = 0 


i.e., the second derivative at the wall is independent of 
the pressure distribution. This suggests the assump- 
tion that there is a universal spanwise layer profile 
0/V, which will represent any v-profile if it is reduced 
to the proper scale. (From the equation of motion of 
the spanwise layer, it also follows that 0*v/0z* vanishes 
at the wall, which gives further support to the above 
approximation.) The boundary condition (3.4) is a 
special case of the condition (1.3) for U’ = 0 orm = 0, 
so that the universal profile should be the special case 
of the family (1.5) for m = 0 (i.e., the Blasius profile), 
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v/V = f(z/0,,m = 0) | 


say, where 


0, = | : (1 — °) de 3.6) 
Do 7 


This is the basic simplifying assumption made by Wild, 
who used the Pohlhausen quartic family of profiles. 
The slope of v at the wall is 


(1 V)(Ov Oz), o> — 6, 


/ 


recalling the definition of / from Eq. (1.5a From 


Eq. (1.11), together with Eq. (1.13), one obtains 


ke o = (1 | = @a/2 = 0.225 


so that the momentum equation [Eq. (3.3)| may be 


written as 
0,(d/dx)(U6,,) = av/2 (3.7) 


The ‘‘mixed’’ momentum thickness @,, [Eq. (5.2)] is 


Put 
6,,/0, = 5S, 0,/0, = q (3.9) 
so that, with the new variable z = @,-¢, Eq. (5.8) be 
comes 
s = | f(¢, m)[1 — g (¢/q) lde (3.10) 
0 
or 
s = s(q, m) (3.11) 


In other words, s is a universal function of g and m; 
the relationship may be evaluated when the family of 


profiles /(¢, m) is specified; and g(¢) = /(¢, 0). Writ- 
ing for the momentum equation [Eq. (3.7 
g,(d/dx)(U0,s) = av/2 (3.12) 


U’, 6,, and m are given functions of x; the latter two 
quantities follow from the calculation of the chordwise 
layer, independent of spanwise flow. As g is a pre- 
scribed function of s and m, Eq. (3.12) is a differential 
or the “‘scale’’ to 
This is 


equation for s (or q)—1.e., for 6, 
which the universal v-profile must be reduced. 
the only remaining unknown after making the basic 
simplifying assumption. 

The function s(g, m) has been evaluated using the 
Pohlhausen-quartic for /(¢, m) in Eq. (3.10), which is 
a straightforward procedure. The results are plotted 
in Fig. 2; the curves drawn there have been inter- 
polated for equal steps of 0.05 of the parameter m. 
This plotting of s against g for different m not only 
greatly reduces the amount of numerical calculations 
involved for the determination of the spanwise layer 
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or, with Eqs. (1.6), (1.10), and (1.18), 
val 
_— a —_ — a , ° . 
gs = a/(L — 2m) aK/m(K + 1) (3,13) que 
Tet ~ ‘ ( for 
Using Eq. (1.13) or (1.19), 
+ eooKxe (Exact) r ture 
_—_-—O— — WIL O-POonLMaUGEN = ” 
ron coonc’s ease 9 ~<~*%°///) qs = a - (6—-1)K +1 (3.14 bec 
ee j a+ (b — 2)m K+ 1 cou 
kno 
Finally, to find the value of g and s, Eq. (3.14) is used, of f 
together with the function s(qg, m), plotted in Fig, 2. T 
The values are given by the intersection points of the sibl 
hyperbolas, Eq. (3.14), in the qs-plane with the 5 any 
curve for corresponding m. The solutions lie on the wri 
chain-dotted line in Fig. 2. The exact values of g and 
s, from Cooke, with the corresponding exact m, are 
also given in Fig. 2; the agreement in the pressure. 
oe drop domain is good but is considerably worse for pres. " 
. q 
sure rise. : 
cori 
Comparison with Sears’ results 
folic 
Another set of values, based on exact calculations, fron 
was taken from Sears’ calculations,” viz., her 4 
m q Ss com 
— 0.0766 1.377 1.970 the 
—0.040 1.303 1.693 \ 
+0 .0922 1.456 1.487 ; 
. . : be ¢ 
Again the values for negative m fit well to the curves of att 
Fig. 2; only the last set is not well represented by the pres 
function s(q, m) of Fig. 2. Int 
= g . ° ent 
g 6 Procedure for Practical Calculation 
Sup 
Fic.2. The function s(g, m) as evaluated from the Pohlhausen It was planned at first to improve the basic functions This 
quartic (solid lines). The points + are evaluated from Cooke's ag. m) of Fig. 2 bv use of exact solutions, but this wa ; 
exact calculation for the velocity distribution U’ = cx*; the HG, We) Os Fig. 4 od wate cunet t aunts a ee eee sure 
corresponding exact values of m are also indicated. The approx- abandoned in view of the scarcity of known exact not 
imate solution for Cooke’s cases, using the Pohlhausen quartic, E 
lie’on the dash-dot line. ;' 6) 
| port 
. nn i ‘ fore, 
but introduces the possibility of a general discussion, 
: ‘ | soug 
as will be seen subsequently. 1-0 oo +508 —— 6 
‘ me repr 
‘ F 9 9 dom: 
Comparison with Cooke’s results 0 
° the § 
First, the method will be checked for the velocity o-6 < ably 
distributions U = cx*. Recently, Cooke" has carried eo i funct 
out calculations, determining the v-profiles in these a? ain k Fig. 
cases by exact calculations, based on the Hartree pro- O6 re eee 
files. His results have been plotted in Fig. 3 against < " LL. O47 
2/0,. If Wild’s hypothesis is justified, the profiles, yg o a Upor 
reduced to this scale, should lie on the same curve for 04 a ‘ “ee multi 
all K. Fig. 2 shows that this is well fulfilled, and the P " 0.43 equa’ 
Blasius profile, which is included in Cooke’s profile- » - Lo 
family for K = 0, may be accurately taken as repre- 2 [- , i ai 
sentative of this universal relationship. - < oe 
7 a as , ‘i —— 
From the similarity of the boundary layers for all x ‘ 0 1 
which exists in the Falkner-Skan-Cooke case, it follows 1 xc 
k , a , 0 2 7. 8 0 2 Uy. 
that s and qg are constants for a given K. Then, from A . 
inn S : 
Eq. (3.12), ¥ 4m 
7 > 3.8 
: et F ; Fic. 3. The velocity profiles of spanwise flow calculated by 6 GF 
qs0,(d dx)(l 0,) = gs(l 6,7 + l 0,0,°) = av/2 Cooke for U = cx*, plotted against 2/0, 7 S$.5 
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The curves calculated from the Pohlhausen- 


values. 
quartic are satisfactory in the domain of pressure drop; 


for pressure rise, on the other hand, it may be conjec- 
tured that the “‘universality’’ of the function s(qg, m) 
becomes questionable anyhow, so that improvements 
could be achieved only when far more exact results are 
known, and they would apply only for a limited “‘type”’ 


of pressure distributions. 

The knowledge of the function s(g, m) makes it pos- 
sible to solve the momentum equation |Eq. (3.12)] for 
any U-distribution by a step-by-step method; if it is 


written in the form 


0.225p « dx — 
5 = - (3.15) 
Ue; Jy gee 


q may be chosen for every interval of integration to 
correspond to the actual value of m(x) and s. The 
value of g for the stagnation point of the chordwise flow 
follows from m(0O) = —0.075. With this value, gs = 3 
from Eq. (3.14). Thus, one can always proceed step 
by step from the stagnation point of the chordwise flow 
component, determining s and thereby the details of 
the boundary layer of spanwise flow. 

Actually, this step-by-step procedure can frequently 
be eliminated in favor of a simpler method. By a 
further simplification, which is satisfactory only for 
pressure drop, a direct analytic solution may be found. 
In this domain of negative m, the curves s(q) for differ- 
ent m lie near together. Actually, the dependence of 
s upon m is much smaller than its dependence upon g. 
This may be seen also from direct reasoning: For pres- 
sure drop, the members of the profile family (1.5) do 
not differ greatly in shape, and, if s is evaluated from 
Eq. (3.10), the ‘‘scale’’ influence—i.e., g—is more im- 
portant than the “‘shape’’ difference—i.e., m. There- 
fore, an analytic dependence of s upon qg alone was 
sought, satisfying the following conditions: It should 
represent closely the actual dependence, Fig. 2, in the 
domain near the stagnation point; and it should follow 
the general trend of the curves for negative m, prefer- 
ably along the line of Cooke’s solutions. A simple 
function with these properties is, as may be seen from 
Fig. 2, 


s=@ (3.16) 


Upon introducing g from here into Eq. (3.12) after 
multiplying both sides by (U/@,)”’, the following 


equation is obtained: 
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$ (U6,s) = 


dx 2 


(U60,s)"* 
or 


2d : 
] (U0.s)"* = 


2 i, av =) 
3 dx 2 WW, 
which may be integrated as 


5s = 


o~ ; 
(U6,) 


3 “Ww + 
6, 4 0 6, 


Thus, Eq. (3.17) yields s by a simple quadrature. 
Within the limits of the additional simplification, 
Eq. (3.16), g is also given: At this point, 
a somewhat different procedure is proposed which 
permits one, by a simple iteration, to avoid almost 
entirely the consequences of this restriction |[Eq. 
(3.16)]. From s, as given by Eq. (3.17), and from m, 
g may be determined for all x from the curves of Fig. 
2. With these new values of g, s may be recalculated 
by Eg. (3.15), yielding a higher approximation. Ac- 
tually, a significant difference between the new values 
of s and those of Eq. (3.17) may be expected only in the 
here, if necessary, the method is 


g= Vs. 


pressure-rise domain; 
easily changed to the step-by-step procedure described 
above. However, Eq. (3.17) is the solution of Eq. 
(3.15) with the approximation g = Vs; even if the 
resulting values of g show some difference when com- 
pared with this assumption, the actual values of s will 
show only insignificant differences from the result of 
Eq. (3.17), s being the result of an integration where 
the integrand contains q. 

Finally, an effort was made to improve the assump- 
tion (3.16) by including the dependence upon m. It 
may be shown that every relationship 


gq = s"F(m) (3.18) 


with an arbitrary constant m and an arbitrary function 
F(m), permits an analytic integration of the momen- 
tum equation [Eq. (3.12)]. The result is 


, 1 x ee 
(U6,s)"t! = m+ i f - dx (3.19) 
2 . F(m) 


Unfortunately, it is not possible to fit Eq. (3.18) with 
enough accuracy to avoid the use of the actual curves 
of Fig. 2 entirely; thus, the iteration method described 
above remains preferable. 


TABLE 1 

Calculation of the Spanwise Boundary Layer for a Circular Cylinder 
1 x,cm 0 1 2 3 4 5 5.5 6 6.5 6.75 6.94 
2 U, cm/s 0 7.106 13.93 20.2 25 29.1 30.2 30.6 30.3 29.8 29.3 
3 O., mm 0.1034 0.1042 0.1062 0.1108 0.1189 0.1338 0.1441 0.163 0.194 0. 2038 0.2224 
4m —0.075 0.0748 —0.0744 —0.0714 —0.0645 —0.049 0.032 —~0.004 +0.037  +0.101 +0. 157 
® S, Eq. (3.17) 2.06 2.05 2.04 2.01 1.96 1.865 1.79 1. 667 1.512 1.495 1.407 
6 4, Fig. 2 1.46 1.46 1.45 1.44 1.425 1.39 1.37 1.34 1.30 1.39 1.45 

2.04 2.02 1.98 1.93 1.83 1.76 1.65 1.48 1.50 1.42 


7 S$, Eq. (3.15) 2.06 2 

















560 JO URNA L OF THE AERONAU TICA L 
4s” 
+7770 17-77 Un GENES 
STAGNATION 4 | er 
LINE | Se 
Y — 
aed . / 
A of 
/ 
| yo | 
| / 














SEPARATION 


A. separation | 


Fic. 4. The stream line near the wall, calculated by the 
approximate method of the present paper (solid line), and the 
potential stream line outside the boundary layer (dotted line) 
on a circular cylinder yawed 45°, plotted on the developed 
cylinder surface from “‘stagnation’’ line to laminar separation. 








Example: Circular Cylinder 


The following example will illustrate the procedure. 
The spanwise boundary layer was calculated for a 
circular cylinder, based on the chordwise velocity dis- 


tribution observed by Hiemenz.? The chordwise layer 


xe. 





Fic. 5. Kerosene-lampblack pattern on the developed surface 
of a circular cylinder, which was exposed to the airstream at 45 
yaw. The calculated stream line near the surface and the 
“stagnation”’ line (top) are plotted. The separation line shows 
clearly on the pattern, as well as a subsequent small region of 
Jaminar flow which is purely spanwise near the cylinder. 


SCIENCES—AUGUST, 1952 
calculations were taken from Pohlhausen’s original 
treatment;*’ the quantities needed for the spanwig 
layer calculation are listed in columns 1—4 of Table | 
Column 5 gives s from Eq. (3.17); the values of qin 
column 6 are found from Fig. 2 with the previous values 
of sand m. The difference of these g-values from the 
quantities V's, calculated from column 5, is not iNsig. 
nificant, but the recalculation of s from Eq. (3.15) shows 
negligible differences with column 5. Here, the pres 
sure-rise domain was small, so that practically no 
further corrections became necessary. 

The slope of the stream line nearest to the wall js 
given by 
dy Ty 


(a/2)(V/6,) _ 0.225 V 
L(U/0,.) ql; U 


= (3.20 
dx fz 


For the stagnation point, the numerical value is 
(1,/Tr)x<9 = 0.46V/U 


that is, the slope at the wall is considerably less than 
outside the boundary layer (where the slope is V/U), 
as was previously found by Prandtl’ and Sears.” 
Near the stagnation point the stream lines at the wall 
and also outside the layer are logarithmic curves, as V 
= const. and U ~ x; the wall stream line lags behind, 
and this behavior may be expected also further down- 
Then, 


approaching laminar separation, the slope at the wall 


stream up to the minimum pressure location. 


increases sharply as /, > 0. 

For the present example, the wall stream line and 
the potential stream line outside the boundary layer 
are drawn in Fig. 4 on the developed surface of the 
The 


evident in these curves. 


cylinder. characteristics described above are 

These theoretical curves have been compared with 
experimental surface stream-line patterns obtained by 
a kerosene-lampblack technique in a wind tunnel at 
Cornell University, and a satisfactory agreement has 
been found, as shown in Fig. 5. The Reynolds Num- 
ber, based on the velocity component normal to the 
cylinder, was definitely subcritical (Re ~ 20,000). 

The question of whether the flow after the laminar 
separation point, and especially the transition mech- 
anism, is also independent of the spanwise flow for in- 
finite cylinders is not quite settled, either in theory or 
by experiment. Loftin® have made 
pressure distribution measurements about yawed ctr 
cular cylinders in the critical Reynolds Number range. 
If the independence principle applies to transition, the 
critical Re Number, based on the velocity component 
normal to the cylinder, should be independent of yaw. 
This rule was found to be confirmed, not exactly but asa 
“rule-of-thumb,” up to 45° Their _ results 
at 60° yaw show a marked difference and reduction 
of Re,,;,.- 


Bursnall and 
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BOUNDARY LAYER IN COMPRESSIBLE 


FLOW 


(4) LAMINAR 


Recently, Ilingworth*!' and Stewartson** have given 
a transformation (generalizing a transformation of 
Howarth”*) by which, under certain assumptions con- 
cerning heat conduction and gas properties, the prob- 
lem of the two-dimensional compressible laminar 
boundary layer for arbitrary pressure distribution may 
be reduced completely to the incompressible case. The 
necessary assumptions are: (a) The surface is heat- 
insulated (no heat transfer); (b) The Prandtl Number 
of the gas is unity; and (c) the viscosity yu varies pro- 
portionally to the absolute temperature 7. Although 
it is felt that for many practical applications these as- 
sumptions may not be adequate, only this case will be 
treated here, with the purpose of combining the Illing- 
transformation with the results 


If any of the assumptions made is 


worth-Stewartson 
hitherto obtained. 
dropped, exact solutions are available only for constant 
pressure, and the approximate solutions for arbitrary 
pressure become much more complicated. 

It is supposed that the isentropic solution outside of 
the boundary layer is known—1.e., U(x)—and also the 
Mach Number .\/(x) is given along the surface. Let 
T) be the stagnation temperature, 7\(x) the tempera- 
ture outside the boundary layer, and 7(x, z) the value 


inside the layer. The isentropic value 7) is given by 


To/T; = 1 + [(y — 1)/2]M? (4.1) 


Inside the layer, the temperature 7’ under assumptions 
is given by the law of Busemann and 


(a) and (b 
Crocco :** 


(4.2) 


On the surface, the stagnation temperature 7) will be 
found for all x, and the viscosity will have its corre- 
sponding value wo. Across the layer, the pressure is 
constant; for a given x, the density p varies according 


to the formula 
pl = pl) (4.5) 


The properties of the incompressible fluid, which 
will be treated after transformation, are normalized, 
for convenience, in the following way: Density and 
viscosity have constantly the values that are found at 
the stagnation point of the compressible fluid—i.e., 
Pi = poy My = pwo—and thus », = vw. The following 
formulas are written out for this case. 

According to the Illingworth-Stewartson transforma- 
tion, which is given here without proof, new coordi- 
nates x;, 2; are introduced, in which the flow problem re- 
duces to the incompressible case, 


* x *\ 3By—1)/2(y7-D 
x; = / (7) dx 


(4.4) 
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8 (TAP 2O-9 78 
oo =) dz = 
To o Pl 
(4.5) 


To a 


In addition (and in contrast to the Mangler transforma- 
tion), here the velocity distribution outside the bound- 
ary layer must also be transformed, by putting 


Ui(x;) _ (To 7)’ . U(x) (4.6) 


at a point x; obtained from Eq. (4.4); actually, U; is 


proportional to M/(x). The same relation holds inside 


the boundary layer between the values u, and u: 
eh . I/e — 
U(X, 21) = (10/71) “u(x, 2) (4.7) 


Consistently with Eqs. (4.4) and (4.6), it follows 
[Using Eq. (4.1)] that 


dU, _ ay I/a-YqGU 
dx; \T; dx 


For the transformed problem, the solution may be 
then, the results must 


(4.8) 


found according to Section (1); 
be transformed back. The following formulas will give 
the results explicitly. 

The momentum thickness in the compressible layer 


has to be redefined as 


| (: _- *.) pu dz 
0 U pil’ 


But the calculations for the transformed incompressible 


0, = (4.9) 


case will furnish the quantity 


0 l i U, 


and, according to Eqs. (4.6), (4.7), and (4.5), these two 
momentum thicknesses are related as follows: 


.~ (7,/to"?"~-" & (4.10) 


The quantity @,, is given from Eq. (1.1) by 


"x 
6,,2 = an,” | U?-' dx 
ih) 


If, here, dx;, U,; and @,, are replaced by Eqs. (4.4), (4.6), 
and (4.10), respectively, the result is 


9 


‘y [v¥+D/Q 1)] b/2) 
o- = avy (3) U 0 «x 
1 
ax 4% 2) 27-1 y¥-1)] 
LG: 


or, for y = 1.4, 


° 3 x 7 1,5 
0,2 = 0.45 (7) U-* (*) Usdx (4.1 1a) 
T\ 0 \To 


This formula gives 6, in the compressible case directly, 
without the necessity of carrying out the transforma- 


tion. 
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For the further investigation of the boundary layer, Note that wo7)/7») = ui, where yu; is the actual viscosity 
it is necessary to determine the shape parameter in the outside the layer, for the assumption hp ~ 7. , 
incompressible case, m;, which makes it possible to use The momentum equation for the compressible lam. “2 
the very relations H;(m;) and 1/,(m;), discussed in  inar boundary layer—viz., - 
a , . — qt 
Section (1). According to definition, de, . ; . = - 
- =e U2— + (0 + 2 — M’)@, l = (4.16 
m, = —(dU;/dx;) (6,;7/ vo) dx dx py & — ‘ 
or, using Eqs. (4.8) and (4.10), (y = 1.4), may be found by transforming Eq. (1.2) by the formy- “ 
dU 0,2 (T,\\8-20/0- las of this section. 
a = = in (7) = To summarize for the calculation of a compressible 
viele . Pe ee laminar boundary layer, the following steps are neces. = 
aU 6," (7) (4.12) sary under the assumptions made: (a) Calculate n 
dx vw \To 6, as function of x, Eq. (4.11); (b) calculate m; as func- 
Now, H, and |; can be determined from the universal tion of x, Eq. (4.12), and read /1/,(m,) and /;(m,); and 
relationships (1.7) and (1.8) for incompressible flow (¢) determine H/ and 7 according to Eqs. (4.14) and 
(Fig. 1, Eqs. (1.11) and (1.13)], and the remaining (4-15), respectively. The reduction to the incompres- wi 
question is the determination of H/ and of the shearing sible case is exact; the momentum method used for val 
stress 7 at the wall in the compressible case. the solution of the latter is naturally subject to restric- 1S 
The displacement thickness for compressible bound- tions discussed in Section (1). 
ary laver must be redefined as It may be mentioned that the Mangler transforma- ¢ 
ae z tion discussed in Section (2) also holds for compres- as 
5* = F ; (: = *) dz = sible boundary layer, so that it may be combined with (5 
0 pil the Illingworth-Stewartson transformation. It is easy, m. 
[ ° (= re *) P ds (4.13) correspondingly, to combine Eqs. (2.6) and (4.11) so te 
4 7: OF te " = that 6, for bodies of revolution in compressible flow may 
; ; , be obtained in one step. The result is se 
where the second form is found by use of Eq. (4.3). su 
Replacing 7/7, by Eq. (4.2), after some rearrangement 9 > _ 9.45, (7) r-2U-8 f GF) or dx (417 he 
it follows that T; ia To at 
aa To - (: e *) P de + Eqs. (4.12), (4.14), and (4.15) remain unchanged. . 
7; te U/ pr For yawed cylinders, on the other hand, the “‘inde- 
T» ‘ u\ up pendence principle’ for the chordwise flow used in Sec- P 
¢ = ) [ (1 - ) 42 tion (3) does not hold (with the exception of iso- Ps 
T . U/ Up thermal fluids), as was discussed by Struminsky.* ” 
The integral in the last term is 6,, Eq. (4.9). Upon Compressibility effects in this case will be subject to me 
dividing by @,, the first term on the right-hand side will further investigations. » 
be recognized as (7)/7\)H;, if Eq. (4.5) is used, so that, Finally, it should be mentioned that an analogy exists st 
finally, one obtains between the problem of the yawed cylinder and the , 
H = (T:/T)H. + (Th/T) — 1 (4.14) heat-transfer problem at low Mach Numbers, discussed : 
by Goland.” Results of Section (3) may be used for 
The wall shearing stress in the compressible flow is this problem. Actually, Cooke’s velocity _ profiles = 
r = po(0u/08),a0 correspond to temperature profiles that, in the thermal ; 
formulation, have been computed earlier by Tani,’ 
where yo is the value for stagnation temperature, the Tjfford,*4 and Oldroyd.” a 
actual temperature of the wall. By use of Eqs. (4.7) - 
and (4.5), 7 may be also written as (5) REMARKS ON TURBULENT BOUNDARY-LAYER \s 
as 1)/(y-1) a CALCULATIONS . 
= fs — : 
G = =O An approximate calculation of the momentum thick- te 
. ; : j ; ness in the case of the turbulent boundary layer (for it 
On the other hand, using the quantity /; introduced in j as 
Shige 5 = i : two-dimensional incompressible flow) may be simply a 
Section (1), Eqs. (1.5a) and (1.16), 7 becomes é ‘ ‘ : : 
derived from the momentum equation, now written ti 
TAO ens OD, in the form 0 
’ vi d6, 4+ (H+ 9) 92 dl Pe. (5.1) tl 
Here, U; and 6,,; can be expressed by Eqs. (4.6) and dx U dx pU? “ 
(4.10), yielding, finally, if the turbulent wall shearing stress 7 is expressed by tl 
t = 1ygio(T1/To)(U/0,) (4.15) the empirical formula* t 
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+/pU? = 0.0128Reo,*** = 0.0128(U@,/v) > (5.2) 


and an “average” constant value of J/7 is assumed; the 
proper empirical value is 7 = 1.4. The momentum 
equation becomes, after multiplication by 1.25U- 


Reo,” 
y-25(1,25U* 58,0290," + 4.25U%%U'0,1-%) = 0.016U! 


or 


(d/dx)(y-**>U**9,!-*) = 0.016U* 


This gives, after integration, 


6, = 0.0366»? U8 (/ U'* dx + c) ~ (5.3) 


where Xo is the value of x at the transition point and C 
must be determined by 6, = 6@,, for x = Xo, where @:, 
is the momentum thickness at the transition point. 

Eq. (5.3), which was first derived by Buri,” is a 
“counterpart” of Eq. (1.1). Although Buri’s methods, 
as a whole, are now generally considered obsolete, Eq. 
(5.3) still can be taken as a sufficiently good approxi- 
mation if only the momentum thickness has to be de- 
termined. 

New questions will arise if the point of turbulent 
separation must also be found. Buri, originally, as- 
sumed a one-to-one correspondence between the mem- 
bers of a turbulent one-parameter family of profiles 
and a shape parameter (which he formed in analogy to 
the laminar case)——namely, (0,U’/U)Re,”* It can 
be computed immediately if 6, is known. Later ex- 
periments have confirmed the assumption of a one- 
paraineter family but not the unique correspondence 
to Buri’s parameter. A new idea was introduced in 
Gruschwitz’s method,?® who gave a purely empirical 
equation for the development of a geometrically defined 
shape parameter. In the newer versions of von Doen- 
hoff and Tetervin® and of Garner,*! the geometrical 
shape parameter is directly 7. Accordingly, the de- 
termination of the turbulent boundary-layer develop- 
ment consists of a simultaneous solution of the momen- 
tum equation [Eq. (5.1)] and the empirical equation for 
the development of 7. Approximately, the two equa- 
tions may be solved separately, taking an average H 
in the momentum equation or the solution to Eq. 
(5.3); then, when @, is given, the empirical equation 
for H may be integrated independently. 

Recently, Tetervin and Lin** have made an attempt 
to interpret the empirical equation as an additional 
integral relationship for the boundary layer—namely, 
as an energy equation or moment-of-momentum equa- 
A discussion of this point of view was carried 
out by Granville.** 

It is not intended here to go into further details of 
the problem in the turbulent case but merely to point 
out that the situation is, in reality, similar to that of 
the laminar case, although the emphasis from the prac- 
For the laminar 


tion. 


tical point of view may be different. 
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problem also, the method based on the momentum equa- 
tion as outlined in Section (1) gave a fairly good formula 
for the momentum thickness—namely, Eq. (1.1). In 
the domain of pressure rise, for the determination of 
the laminar separation point, the one-to-one corre- 
spondence of the shape parameter m and of the particu- 
lar member of the profile family which gives best repre- 
sentation to the actual profile becomes questionable. 

Attempts have also been made to obtain improve- 
ment in the laminar case by the use of the energy equa- 
tion. While in the turbulent case even the reinterpre- 
tation of the additional empirical relationship accord- 
ing to Tetervin and Lin still must contain empirical 
functions, the energy equation for laminar flow gives 
an additional relationship without empirical assump- 
tions, because the boundary layer profile and the shear- 
ing stress are connected by Newton’s law. First 
Wieghardt** has used the energy equation, in addition 
to the former assumptions, and has based his method 
on a two-parameter family of profiles. The improve- 
ment was obtained at the cost of greater complication. 
Walz®*® has proposed to use a one-parameter family 
again but to drop the fulfillment of the boundary con- 
dition (1.3) and thus to abandon the one-to-one corre- 
spondence of the actual profile with the parameter m 
and to use the energy equation instead. This caltula- 
tion is, in principle, analogous to the procedure for tur- 
bulent layers, except that no empirical assumptions on 
the shearing stress are necessary. For Walz’s method, 
much better agreement with exact solutions and with 
experiment is claimed when the laminar layer develop- 
ment is investigated for pressure rise. The improve- 
ment is explained by the fact that an integral relation- 
ship fixes a member of the profile family in a much 
more insensitive way—i.e., much more independently 
of the actual and accidental properties of the profile 
family—than does the parameter m, which is connected 
with the second derivative at the wall, Eq. (1.5a). 

The conclusion may be drawn that the approximate 
boundary-layer calculations based on the momentum 
equation alone give simple and useful formulas for the 
momentum thickness, both in laminar and in turbulent 
flow; but for the question of separation, an improve- 
ment may be necessary by the use of two integral rela- 
tionships as discussed above. The practical difference 
between the laminar and the turbulent case is that the 
question of laminar separation is, in reality, over- 
shadowed by the transition problem (Schlichting*), 
while the development of turbulent layers in domain of 
pressure rise and their separation are of the greatest 
practical importance. 

For the turbulent boundary layer on a body of revo- 
lution, starting from the modified momentum equation 
[Eq. (2.1)], the counterpart of Eq. (2.6) was established 
by Truckenbrodt.** With the same assumptions that 
led to Eq. (5.3) and by analogous steps, only taking 
account of the difference in the momentum equation, the 
following formula may be proved: 
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6, = 0.0366y—°-?U—*r—! Cf Ur! dx + c)" (5.4) 


Recently, it was found by Young and Booth,* sup- 
ported by experimental evidence, that the ‘“‘independ- 
ence principle’ for yawed infinite cylinders, as stated 
in Section (3), is also valid for the turbulent boundary 
layer. According to this principle, the chordwise layer 
may be determined as if the spanwise flow were not 
present. Altman and Hayter* gave further experi- 
mental evidence to this result for the turbulent boundary 
layer on a swept wing. 

The statement of Young and Booth is, in some re- 
spects, surprising and is by no means trivially identical 
with the principle in the laminar case. To clarify this 
point, a simple example will be discussed—namely, the 
case of zero pressure gradient on the yawed flat plate. 

In this case, it was found for the laminar boundary 
layer by Prandtl,’ Jones,'® and Sears” that the veloc- 
ity profiles have similar chordwise and spanwise com- 
ponents, both being similar to the resultant profile, and 
that all velocities are parallel to the resultant stream- 
Thus, for the momentun 
= §,, = 0 or that 0 
Consequently, 


direction outside the layer. 
thicknesses, it follows that 6, = 0, 
is invariant with respect to direction. 
for zero pressure gradient, it is easily seen that the lam- 
inar boundary layer may be treated in either of two 
equivalent ways: (1) Follow the resultant stream di- 
rection and disregard the yaw of the leading edge: and 
(2) follow the chordwise direction and disregard the 
spanwise flow. Only the second way uses the inde- 
pendence principle, but the two methods are equivalent 
in this special case. 

For the turbulent boundary layer, these two points 
of view cannot be equivalent. Assume again that the 
resultant velocities through the boundary layer are 
parallel to the main stream direction (as the experiments 
of Young and Booth indicate). If the first method 
(disregard of yaw) is accepted, then the resultant wall 
shearing stress, computed from Eq. (5.2), is 


onan. 0% /*%/ TT? 79) 1.75 ‘a 
Tres, = 0.0128pv—* 9" *(V U? + V2)" = (5.5) 
and the component in the x-direction is evidently 


r, = 0.0128pv—°-259-%U(V U2 + V2) (5.6) 
On the other hand, if the spanwise flow is completely 
disregarded, as the independence principle requires, 
and Eq. (5.2) is used to give the x-component of the 
wall shearing stress directly, the value of 7, is given by 

7, = 0.0128 py—%- 2590-2 [1-7 (5.7) 
The two results differ; this is evidently the conse 
quence of the fact that the turbulent shearing stress is 
proportional to the 1.75th power of the velocity, so that 
it makes a difference whether a component is taken of 
the resulting shearing stress or a shearing stress com- 
ponent is calculated from a velocity component. In 
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the laminar case, where 7 is proportional to the first 
power of the velocity, this difference does not exist. 
At first instant it seems to be more natural to accept 
the former point of view, Eq. (5.6), but experiments of 
Young and Booth® are definitely in favor of the second 
), and thus in favor of the in 


point of view, or Eq. (5.7 

dependence principle. The resultant shearing stress 
calculated back from Eq. (5.7), becomes 

= 0.0128py—0->9-V U2 + V2U% (5.8 


Tres 


In their paper, Young and Booth*® have simply 
stated the ‘‘evident”’ independence of 7, from the span 
It is felt, however, 


wise flow, or the use of Eq. (5.7). 
that this point needs further clarification from the 
theoretical point of view. For laminar flow, the inde 
pendence principle was an exact conclusion from the 
boundary-layer equations; only the assumption of an 
infinite cylinder represents a certain idealization. In 
the turbulent case, on the other hand, the Reynolds 
shear stresses are not a priori connected with the mean 
flow quantities, and the statement of the independence 
principle implicitly contains a hypothesis on the turbu- 
lent transfer mechanism. It is believed that the fol- 
lowing picture may help for the visualization of the re- 
sult of Young and Booth: The (resultant) growth di- 
rection of the turbulent boundary layer is, for geo- 
metrical reasons, the chordwise direction and not the di- 
rection of the resultant flow. It can well be imagined 
that, for the establishment of a mechanism of turbulent 
momentum transfer, only that velocity component plays 
a role which lies in the direction of boundary-layer 
growth—1.e., U. The spanwise component V is paral- 
lel to lines of constant boundary-layer thickness, so 
that it appears possible that the turbulent transfer 
mechanism becomes independent of this quantity 
This may be helpful in understanding Eq. (5.7). 

From this point of view, it seems to be nontrivial 
that the boundary-layer profiles (for zero pressure 
gradient) come out similar in the U- and V-direction 
On the other hand, it may be recalled that, once a 
mechanism of turbulent transfer is established, it will 
provide the transfer of any quantity—e.g., heat or 
particles—in the same way as the transfer of momen 
tum. In the present case, the transfer of V-momentum 
becomes similar to the transfer in the l/-direction. 

If the “independence principle’ is understood for 
zero pressure gradient, its generalization for any pres 
sure distribution appears natural. For the spanwise 
flow, the momentum equation | Eq. (3.1)] or 


pV(d dx)(U0.,) = Ty (5.9 


holds; 


for r,, it cannot be used in the same sense as in the 


Young and Booth* assume a ‘‘univer- 
then, making an assump- 


laminar case. 
sal’’ spanwise flow profile; 
tion on the thickness of the spanwise layer, they use 
the momentum equation for the determination of 7 
In the laminar case, 7, was related to the spanwise 





but, because no simple hypothesis is available 
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flow profile, and the thickness of the spanwise layer was 


computed. 
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Lateral Buckling of Nonuniform Beams 


Bruno A. Boley and Vincent P. Zimnoch 

Department of Aeronautical Engineering, The Ohio State 
University, Columbus 

April 22, 1952 


T# NUMERICAL METHOD PRESENTED HERE was first developed 
in reference 1 in connection with an investigation of the lateral 
buckling of wind-tunnel load members. Such members are 
usually tapered; for this reason the methods of reference 2, em 
ploying either differential equations or energy considerations, 
were found to be cumbersome. Applications of these to the 
problem of lateral buckling have recently appeared in references 
3 and 4. The approach adopted in reference 1, and described 
below, makes use of the procedures employing influence co 
efficients developed in references 5 and 6. 

The beam in question is replaced by an imaginary beam made 
up of a number of sections, each with uniform properties. Each 
end point of these sections, called a ‘‘joint,’’ has six degrees of free- 
dom—namely, the displacements u, v, and w in the x, y, and s 
directions, respectively, and the rotations a,, a,, and a, about 
the axis indicated in the subscript. Buckling is assumed to take 
place in the direction and under the applied loading 7 shown 
in Fig. 1; under these conditions the quantities w, v, and a, may 
be neglected during the buckling process. Three equilibrium 
equations may then be written at each joint, corresponding to 
the remaining three degrees of freedom. At the ith joint, each 
of these equations is of the form 
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> (a,ju; 4 bijayj + Cijazj) = 0 (1) 


j 1 


where the quantities a;;, b;;, and ¢;; are influence coefficients and n 
is the number of joints in the structure. The motions of the jth 
joint are denoted by u;, ay;, and @,;. A set of (3) simultaneous 
homogeneous linear algebraic equations thus results; from it, 
the buckling load can be calculated by the methods of references 
5 and 6. The criterion used here states that the determinant 
formed by the coefficients of this set of equations vanishes at 
buckling. Evaluation of this determinant (which is symmetric) 
is simplified by the method of reference 7. 

The influence coefficients a;;, bj;, and ¢;; are found from con- 
sideration of a single section of the beam. For any such section 
the governing equations are, in the notation of Fig. 1, 


Bu" = Ma, + FAL — s)a, + M, | 
(2 
Ca,’ = —M(u’ — w’) + M, { 


where B and C are the flexural and torsional rigidities, respec- 
tively, M is an applied moment, and primes indicate differenti 


ation with respect to z. Eqs. (2) may be combined to yield 


a,” + k’a, = —(R?F,/M)(L Z) k™M,/M) (38a) 


where 
k? = M?/(BC) (3b) 


The force F, and the moments /, and j/,, due to a motion of a 
joint in each degree of freedom separately, are found from the 
solution of Eq. (3a) under the appropriate boundary conditions. 


The results are presented in Table | (note that u’ = a, 





TABLE | 
a ay u 
Vi C/L —M U/L 
M, —M b a 
F M/L a ( 
where 
k?B sin kL 
- kL(1 + cos/RL) — 2 sin kL 
— kB sin kL (kL cos RL — sin kL) 
~ (1 — sin RL) [RL(1 + cos RL) — 2 sin kL] 
P — 2k? B sin kL 


~ RL? (1 + cos kL) —2L sin kL 


The method described above was used in a number of examples. 
In each case, Eqs. (1) were set up with the aid of the influence 
coefficient of Table 1, and the determinant of the system was 
equated to zero. The buckling loads obtained for the various 
cases are collected in Table 2. It may be noticed that satisfac- 
tory results are obtained in all examples when six degrees of free- 
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TABLE 2 
w (2) (ED) () (5) 
STRUCTURE BUCKLING LOAD BUCKLING LOAD TOTAL NUMBER | PERCENT ERROR 
AND LOAD BY PRESENT METHOD | BY METHODS OF OF DEGREES | (3)- x 100 
OF FREEDOM Ey) 
UNIFORM ; — 
CANTILEVER, wu » 7 V8 .mVec 3 ° 
MOMENT MAT ee. FT 
END 
e. o TYRE 3 207 
u van 
CANTILEVER, p » #013Vec a 
CONCENTRATED Vas er L 
LOAD P AT END e. - 22tec 6 29 
c 
TAPERED PR, ©3340 LB 3 300 
CANTILEVER, 
CONCENTRATED 
LOAD P AT END 
ee VARIATION e +4260 LB 
Pp. +4740 18 6 1 
s.s* (BC ASSUMED CONST.) 
(Bc) = 7.0 x 10° =, 
tp IN se 
i | = 4400 LB. 5 70 
(BC), + G (BC) 34 e 
®, = 2550 LB 3 26.0 
SAME AS ABOVE, saiiins 
EXCEPT 
R,: 
isucdiea P+ 3460 LB 6 67 
(2 ®8C ASSUMED CONST 4 
(Bc), (7%) (@C eet 
a 
®, #3220 1.8 6 7.0 

















dom are used and that all the buckling loads calculated by the 


present procedure are on the safe side. 
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Prestressing Rectangular Thin Plates 
to Increase Their Buckling Loads 


Chi-Teh Wang, Arthur L. Ross, and Edward L. Reiss 


Professor of Aeronautical Engineering, Instructor in Mechanica, 
Engineering, and Graduate Student in Aeronautical Engineering 
respectively, New York University 

March 12, 1952 


wr A FLAT RECTANGULAR PANEL is under compression in 
the plane along two opposite edges it buckles when the com- 
pressive load reaches a certain magnitude. This buckling load 
depends on the dimensions of the plate, the edge-support condi- 
tions, and the mechanical properties of the material from which 
the plate is made. Thus, for a given panel, to increase the buck 
ling load it is necessary to increase the thickness of the plate if 
the same material is to be used. This, however, increases the 
weight of the panel, which may not be desirable. A more effi- 
cient method of increasing the buckling load without increasing 
the weight is obtained by prestressing the plate. 

It can be shown that, if tensile forces are applied in the direc 
tions perpendicular to the compressive forces, the buckling load 
of the same panel can be increased appreciably. For panels such 


as those usually used in aircraft construction, a simple method 
for producing such in-plane tensile forces is as follows: The flat 
plate is first cold-rolled into a curved plate in the direction per- 
pendicular to the loading. The resulting curved plate wil] then 
be attached to the frame, becoming flat again. By such a proc 
ess, in-plane tensile forces perpendicular to the direction of load- 
ing will be produced by the large deflection of the plate The 
resulting flat panel will then have a much higher buckling load 
To verify the above proposal, a preliminary experiment was 
carried out A flat 19'/, by 14 in. 24ST aluminum plate with a 
thickness of 0.020 in. was first cold-rolled into « curved plate 
with an approximate radius of curvature of 9! , i This re 
sulting curved plate was then clamped into a frame ( Fig. | The 
deflection of a similar plate was measured by a dial gage. It was 
found that the maximum deflection was about one-half the thick 
ness of the plate and the panel could therefore be regarded as 


flat. Strain gages were attached on both sides of the plate at 
9 


four stations (Fig. From the strain-gage readings, it was 


found that the strain in the direction of loading was practically 
zero, while the in-plane tensile strains ey at the station 0), 1, 2 3 
were +60, +80, +10, +15 microin. per in., respectively. With 
these in-plane tensile strains, the buckling load can be computed 
by means of finite difference approximations The nondimen 
sional buckling load .V,6?/D in this case was found to be 210 
By a similar procedure, the corresponding buckling load for the 
same flat plate without prestressing was 101, where D is the 
flexural rigidity of the plate and + is the length of the loaded edge 
This indicates that prestressing will increase the buckling load a 


considerable amount. 
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The fact that the buckling load of a curved plate is increased 
materially by first overcurving the plate was first recognized by 
Welter.2 It was suggested,® 4‘ however, that in the case of a 
curved panel the reason for raising the buckling load is due to 
the favorable deflections caused by the process. There is no 
doubt that the favorable deflection was a contributing factor 
toward the increased buckling load. However, it seems that the 
initial stresses induced by such a process must be equally im 
portant factors if not more important ones in explaining this 
phenomenon This last statement is based on the fact that, in 
the limiting case of a flat plate, there cannot be any favorable 
deflections, and any increase of the buckling load must be caused 
by the initial stresses. More experiments are being conducted at 
the Guggenheim School of Aeronautics of New York University, 


ind the results wil! be presented in a later report. 
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End Correction for Flight-Tested Frequency 
Response Obtained by Laplace Transformation 


G. M. Andrew 
Edwards Air Force Flight Test Center 
April 7, 1952 


hem RECENT FLYING TECHNIQUE for obtaining the frequency 
response of the airplane from the flight test reduces (most 
frequently) to the application of control surface motion as a 
step or pulse input. This motion is recorded as a time history. 
The response of the airplane representing linear and angular 
displacements, velocities, and accelerations are also recorded as 
time histories. The most typical records are shown in Figs. 1, 2, 

3, and 4. 
Fig. 2 represents a heavily damped mode where at a certain 
h(T). For Fig. 3 at time 


T, the slope becomes constant (1) = constant 


time 7h(t) assumes constant value 


Fig. 4 represents 
a lightly damped mode, where from time 7, 


h(t) = Are~™ sin (qt + yw) 


It is also possible to have a combination of these three basic cases 
The frequency response H(jw) is obtained from time records 


h(t) by Laplace transformation 


H(s) -{ 
0 


$ = ¢ 


e~“h(t) dt 


+ jw 


lore = QVors = jw 


H( jw) - | "¢ 
0 


This is done by numerical integration of the time record 
If the record does not end with value of A(t) = 0, the end cor- 
That means that the numerical integra- 


Joth(t) dt (la) 


rection must be applied. 
tion is carried out up to time 7, and from then on an analytical 


correction isadded. Thus, Eq. (la) is replaced by 


. 
H(jw) = > e~Ieh(t)At + AH(jw) 


(1b) 
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where 


AH(jw) = J, e—Sety(t) dt 2) 
The real part and imaginary parts of Eq. (2) have to be added to 
Re and Jm solution of the first part of Eq. (1b) 
(1) Step 


From time ¢ = 7 (Figs. 1 and 2), 


h(t) = h(T) = const 


Laplace transform of Eq. (3) is (reference 1, transform No. 1.01) 


h(T)/s 
As the step does not happen at ¢ = 0 but with time lag ¢ = 7, 
the total transform is (reference 1, transform No. 10 
AH\(s) = e ~8? [h(T)/s] 3b) 
Substitution s = jw yields, immediately, 
; wp MT) 
AH\(jw) = e~9* : 
jw 
_.. sin wl + j cos wT : 
= —h(7 } 4) 


WwW 
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(2) Ramp 
From ¢ = T (Fig. 3), h(t) has constant slope 
h(t) = const. (5a) 
The Laplace transformation of Eq. (5a) with time lag T is 
(reference 1, transform Nos. 2.101 and 10) 


AHAs) = e~57 [h(T)/s?] (5b) 
Substitution s = jw gives 
AHAjw) = boil [h(T)/(jw)?] 
= —[h(T)/w?] (cos wT — j sin wT) (6) 
(3) Sine Line 
From time t = 7, 
h(t) = Are~™ sin (qt + ) (7) 


Fig. 4 gives the notations, 47, P, A;, Ax; Atand 7 are measured 
from the record. 


Then, the angular frequency: 
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damping exponent : 
n = (In A; — In Az)/At 9 
iead angle: 
y = 2 (1/P) (10 


Expanding in Eq. (7) the term sin (qt + W) gives 
h(t) =A Te nt (cos yp sin gt + sin y cos qt) (lla 


Laplace transforms of Eq. (11) with time lag T yield (refereng 
1, transform Nos. 1.301, 1.3011, and 10) 


gcos Y + (s + n) sin y 
AH;(s) = Ares! : ° : 
(s + n)? + q? 
Substitution of s = jw gives 
; _; 7 fqceos ¥ + (jw + n) sin y 
AH;( jw) = Are jo ; - ; (11b 
(jo + 2)* + q* 


Substitution of numerical values into (11b) allows fast ration 
alization by means of a slide rule, giving Re and Jm parts of 
AH (jw). 


If desired Re and Im part of Eq. (11b) can be computed from Eqs. (llc) and (11d), which were obtained by proper reworking of Eq 


(11b),? 
: A|nsin (wT + y) + (w + g) cos (wT + y) n sin (wl — WY) + (w — gq) cos (wT y) 
Re [AH;(jw)] = ee eee - . (11¢ 
2 n* +- (w + q)® a” -- (o.— gr 
: Ar|ncos (wl + ~¥) — (w + g) sin (wT + wy) n cos (wT — yw) — (w — g) sin (wT y) 
Im |AHdjw)] = T cos (w Y w +g _ w o. . : E (114 
3 n? + (w + q)? n* + (w — g)* 


Selecting y = O or yy = 2/2 simplifies considerably Eqs. (11c) and (11d). 
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Heat Transfer Near the Forward Stagnation 
Point of a Body of Revolution 


M. Sibulkin 
U.S. Naval Ordnance Laboratory, White Oak, Silver Spring, Md. 
April 11, 1952 


ie ORDER TO DETERMINE the temperature distribution over a 
body moving through the atmosphere, a knowledge of the 
local heat-transfer coefficients is required. For slender sharp- 
nosed bodies, the heat-transfer coefficients are frequently ap- 
proximated by using the comparable flat-plate values. 
ever, for blunt-nosed bodies, flat-plate solutions are not appli- 


How- 
cable near the forward stagnation point. Since the greatest rate 
of heat transfer may occur at the forward stagnation point, its 
In this note a theoretical solution 
is given for the heat transfer near the forward stagnation point of 


value should be investigated. 


a body of revolution assuming laminar, incompressible, low-speed 
flow. The comparable solution for two-dimensional flow has 
been given by Squire.'! In the case of a blunt-nosed body moving 
with supersonic velocity, the flow behind the central portion of 
the bow wave is subsonic, and it is possible that a low-speed solu- 
tion, using as ‘‘free-stream’’ conditions those behind the center 


of the bow wave, will apply near the stagnation point. 





Near the forward stagnation point, the velocity, U’, just out 
side the boundary layer may be put equal to Bx, where x is the 
distance from the stagnation point along the surface and 8 isa 


constant. Then the momentum equation is given by 
Ou Ou x O7u 
u +? = fx + v (1 
Ox oy Oy? 


where y is the distance perpendicular to the surface, u and 2» are 
the velocity components in the x and y directions, and +r is the 
kinematic viscosity. The energy equation, neglecting the pres 
sure and dissipation terms by restricting the analysis to low 
speed, is 

oT 


oT oT 
v = K 
Oy? 


q = 


u + 
Ox oy 


where 7 is the temperature and x is the thermometric conduc 
tivity. The continuity equation for axially symmetric flow is 
given by 

[O(ru)/dx] + r(dv/dy) = O (3 


where 7 is the distance of a point on the surface of the body from 
the axis of symmetry. 

Near the stagnation point of a blunt body, 7 is nearly equal to 
x, and Eq. (3) can be satisfied by putting 


1 dy 1 dy 
“= ; = 
x Oy x Ox 
Now define 
n = (8/v)'*y, W = (Br)'/2x?f(n) 


Performing the required differentiations and substituting in Eq 
(1) gives 

ft — off" =a ey” (4 
This equation has been solved,? and values of f, f’, and f” have 
been tabulated as functions of 7. 
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TABLE | 

o c 0.76300.4 
a 0.6 0.625 0.622 

0.8 0.700 0.698 

1.0 0.763 0.763 

2.0 0.988 1.01 

10.0 1.76 1.92 

Now assume 
T= T, — (T. T « )g(n) 


where the constants 7, and 7. are the wall and free-stream tem- 
peratures Performing the required differentiations and sub- 
stituting in Eq. (2) gives 

g” + 2ofg’ = 0 (5) 
where o is the Prandtl Number and, with the boundary condi- 
tions, 

g=Oatyn=0 


g=latyn= @ 


Integrating Eq. (5) and using condition (6a) gives 


: en a bl 
g= Co) f exp | —2e¢ f dn) dn (7) 
0 0 


Applying condition (6b) gives for C(a) 


' * mn 
: = | exw ( —2¢e / f in) dyn (8) 
Cla) 0 J0 


Eq. (8) was numerically integrated for several values of Prandtl 
Number using the known values of f; the results are given in 
Table 1 2.0, C(@) can be approximated 
by 0.7630" 


The rate of heat transfer is given by 


Nu, k(T Pia (: r) 
q = = —k (9) 
L Oy /y=0 


where g is the rate of heat transfer per unit area, k is the thermal 


In the range 0.6 < oa « 
1 


conductivity, is a representative length, and the Nusselt Num- 
ber based on L, Nuz, is defined by Eq. (9). Differentiating Eq. 
7) and substituting » = 0 gives g’(0) = C, and from Eq. (9) 


1 


Nuy = L(B/v)/*2'(0) = CL(B, v)' (10) 


D, where U' is the free-stream ve 
Setting Z equal to D and 


6 = 3U. 
locity and D is the sphere diameter. 
using the approximate value of C in Eq. (10) gives the heat trans- 


For a sphere, 


fer at the stagnation point of a sphere as 


1 


Vup = 0.763809 (3U2D/v)'/? = 1.32Rp°*o4 (11) 


where Vup and Rp are the Nusselt and Reynolds Numbers based 
upon D. For other blunt-nosed bodies of revolution, the value 
of 8 and its relation to L must be determined from potential 


flow theory or from experimental measurements 
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On Slender Wing-Body Theory 

John R. Spreiter 

Ames Aeronautical Laboratory, N.A.C.A., Moffett Field, Calif. 


April 17, 1959 


een NOTE BY Mires! has called attention to a disparity 
between the results of slender wing-body theory as given by 


myself? and by Ward’ and Stocker.‘ The ideas involved are not 
peculiar to slender wing-body theory but apply equally in the 
complete linear analysis of the same problems. Since Miles’s 
comments are brief, it is felt that a further statement might be in 
order. 

The linear theory of compressible flow involves the solution 
of the Prandtl-Glauert equation for the perturbation velocity 


potential 

(1 — Me)ece + yy + Gu = 0 (1) 

together with appropriate boundary conditions. [Slender wing 

body theory referred to above involves the simplification of Eq 
(1) to the form 

¢yy + ¢2z = O (2) 

but this is of no particular consequence for the following remarks. | 

Upon solving the given boundary-value problem, the pressure 

coefficient Cp for wings (and, in former times, for bodies of revolu 

tion ) is calculated therefrom by means of the linearized Bernoulli 

equation, 
P — Po x 
ais — ( 


Cy, = 2 
y (p0/2) Uo? Uo 


w 


It was later recognized in the theory of axisymmetric bodies that 
certain additional terms should be retained; thus, Lighthill® 
showed that it is consistent with the approximations of the theory 
to let 


Cp = —2(¢2/U0) — [( ey? + 27)/ U0?) (4) 
The reason for the inclusion of the additional terms stems essen 
tially from the fact that ¢, may become considerably smaller 
than gy and ¢,. 

It is now apparent? that a similar situation exists in the treat- 
ment of low aspect ratio wings, as well as slender wing-body 
combinations, and that the pressure coefficient should again be 
computed from Eq. (4). It should be noted that, although the 
velocity components due to thickness and angle of attack can be 
determined separately and superposed to give the final value 
superposition of the pressures is forbidden by the nonlinear char 
acter of the pressure-velocity relation. 

Although the inclusion of the additional terms in the pressure 
relation always changes the pressures, the loading or difference 
in pressure between the upper and lower surfaces may or may 
not be affected. Thus, consider a wing-body combination, and 
let v% and vq represent the v-velocity components due, respec- 
tively, to thickness and angle of attack. The directions and 
magnitudes of 7 and vq are as shown in Fig. 1. The v? term in 
the pressure relation then contributes to the loading as follows: 


vy? — vy? = (% + Ve)? —(% — Va)? = 4a (5) 


where u and / refer, respectively, to conditions on the upper and 
lower surfaces. From such corsiderations, it is apparent that 
the loadings calculated by means of Eqs. (3) and (4) are, in gen- 
eral, different but become the same when the body is cylindrical 
(when 7, = 0). It should be remarked that the determination of 
vy, requires the use of Eq. (1) rather than Eq. (2), as sometimes 
suffices for vg. At the present time, formulas for determining the 
velocity components due to thickness for slender objects of arbi- 
trary shape have been given only for supersonic flow.* ‘ The 
analogous problem for subsonic flow has been treated by Heaslet 
and Lomax in their forthcoming Princeton Series article. 
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The slender wing-body theory of reference 2 was developed 


through the use of Eq. (3) rather than Eq. (4). Since the equa- 
tions are for the loadings and not the pressures, the results are 
nevertheless correct for the wing alone, body alone, and the wing 
mounted on a body that is cylindrical along the wing root; but 
not so for the more general case where the body diameter varies 
along the wing root. 

The specific example of a triangular wing mounted on a conical 
body treated in reference 2 [and independently in reference 6 
using Eqs. (1) and (3)] illustrates the contribution of the square 
terms. The ratio of the lift of the combination to the lift of the 
wing alone L/L,, is plotted as a function of the ratio of the body 
radius to wing semispan in Fig. 2. The results given originally in 
references 2 and 6 are shown by the dotted line. 
shows the values obtained using Eq. (4), together with the req 
uisite additional terms in the potential, to account for the flew 
The solid curve 


The solid line 


around the combination at zero angle of attack. 
was given first by Ward.* The source of the discrepancy was 
not immediately evident because the forces are determined by 
pressure integration methods in references 2 and 5 and by mo- 
mentum considerations in reference 3. The writer wishes to 
thank A. E. Bryson, Jr., and R. H. Edwards, of the Hughes Air- 
craft Company, for their help in clarifying this situation. 

The same considerations also apply to slender cruciform wings 
inclined simultaneously in pitch and yaw. For example, in refer- 
ences 7 and & the loading on the horizontal wing of a cruciform 
arrangement of low aspect ratio triangular wings is found through 


application of Eqs. (2) and (3) to be 


(7) ta|(ds/dx) + B(y/s)] (6 
— = ») 
q /H 


2 


V1 — (y?/s?) 
If the pressures are computed using Eq. (4), the following result 

is obtained 

Ape _ ta(ds/dx) 
q /H 7 V1 — (y?/s?) V1 — (y?/s?) V1+ (t/g?) 


ta By/s l 


where s and ¢ are, respectively, the local semispans of the hori- 
zontal and vertical wings and @ and @ are the angles of attack and 
sideslip. In this case, the additional terms affect only the asym- 
metrical portion of the loadings, thereby leaving the total lift 


and pitching moment unchanged. Although the rolling mo- 
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ments are affected in general, it should be remarked that the total 

. . ~ . al 
rolling moment on a cruciform configuration having identical 
wings is still zero to the accuracy of the theory 
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Comments on ‘‘A Study of Transonic Gas 
Dynamics by the Hydraulic Analogy”’ 


W. C. Randels 
U.S. Naval Ordnance Test Station, Inyokern, China Lake, Calif 
April 24, 1952 


ge tapas | STATES THAT the hydraulic analogy is only valid 


» 


for a water depth of 0.25in. This conclusion does not seem 


to be justified by the discussion. The reasoning depends on the 
distinction between wave velocity and group velocity. How 
ever, in the expression given for the group velocity, c — \(dc/dd), 
the term \(dc/dX) ~ 0 as \ ~ ©&, so that for large \ the group 
velocity is equal to the wave velocity. 

This agrees with the experimental results cited in reference 2 
In this report it was possible, by choosing a large enough model, 
to get agreement between measured and calculated oblique 
shock wave angles over a considerable variation of water depths 
The experimental scatter of these measurements is large so that 


some effects may be masked. 
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Comments on ‘‘Remarks on Dynamic Loads 
in Landing”’ 


E. H. Kramer 
Dynamics Engineer, Wright-Patterson Air Force Base, Ohio 
April 30, 1952 


IX ADAM ZAHORSKI’S PAPER “Remarks on Dynamic Loads 1 
Landing” (JoURNAL OF THE AERONAUTICAL SCIENCES, Vol 
19, No. 4, pp. 258-264, April, 1952), the statement is made 
that other methods ‘‘as proposed by Wright Air Development 
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are simple in calculation but are based on arbitrary as- 


Center, : é ; 
ymptions and seldom agree with experimental results. It is 
onsidered that this statement is misleading, in that Mr. Zahorski 
cons 


did not list the arbitrary assumptions to which he objects, nor 
did he offer any verification that these computations do not always 
agree with experimental results. . 

The basic assumption in the W.A.D.C. method of dynamic 
analysis is the trapezoidal forcing function of which Mr. Zahorski 
apparently approves, and it should be pointed out that, when the 
W.A.D.C. method (A.M.C. Report MCREX4A5-8-2) was de- 
veloped, it was checked against experimental results available 
at that time and was found to give satisfactory correlation. 

I believe Mr. Zahorski’s paper is of value to the structural engi 
neer, for, although the difficulties of dynamic analysis are known 
only too well to the dynamics engineer, most other engineers do 


not appreciate the complexities of this problem 


On Correlation of Vibration and Buckling 
Theory of Plates 


George Sonnemann 
Assistant Professor, Mechanical Engineering Department, Drexel 
Institute of Technology, Philadelphia 


April 25, 1952 


N A RECENT ARTICLE IN THIS JOURNAL, Lurie! pointed out that 
| a correlation between the vibration problem and the buckling 
problem should be feasible. The equations offered in his article, 
though they show an apparent correlation, do not truly represent 
the similarity of the two problems. This must be attributed to 
the fact that a plate will always vibrate with at least a half- 
wave in the .Y-, as well as Y-, coordinate direction. In Lurie’s 
paper the frequency is only a function of the wave length in the 
X-direction. 

In order to establish a more general correlation, it is best to 
confine one’s effort to the energy solution. By means of this 
method it is possible to show that correlation does exist, in gen- 
eral. 

The correlation equation for the simply supported plate is 


wo= Vy D phb* 32k 


where k is the appropriate plate buckling constant. For a simply 


supported plate, the value of & is given by 
k = [(m?/B?) + n?] 


where 8 = a/b. The theory and extension of the correlation for- 


mulas will be reported in the near future. 


REFERENCE 


‘Lurie H. Vibration of Rectangular Plates, Readers’ Forum, Journal of 


the Aeronautical Sciences, Vol. 18, No. 2 p. 139, February, 1951 
+ 


Equations for Estimating Laminar Boundary- 
Layer Effects 


Ralph Ashby Burton 

Teaching Fellow, Department of Mechanical Engineering, 
The University of Texas, Austin, Tex 

May 7, 1959 


I ORDER TO PRODUCE simple integral equations for boundary- 
layer thickness, shear, and separation, the boundary layer 
€quation of Prandtl must be simplified. This equation may be 

written in the following form: 
du 1 dP O7u (1) 


v 9 
dt p dx oy? 
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Fig. 1 


where u is the velocity component in the boundary layer, parallel 
to the surface; x is the coordinate along the surface measured 
in the direction of the free-stream velocity adjacent to the bound- 
ary layer; y is the coordinate measured normal to the surface; 
P is the pressure acting at the surface, and /is time. The quan 
tity p is the mass density of the fluid and » is the kinematic vis 
cosity. 

One method of simplification of the differential equation is to 
linearize the acceleration term in the following manner 
1 dP 


= - p 


pdx Oy" 


V Ou 
2 Ox 


O7u 


is the velocity of the frictionless flow immediately ad- 
While this linearization owes its 


Here, V 
jacent to the boundary layer. 
origin to the precedent of successful similar linearizations in vari- 
ous flow problems, its principal justification is that it permits the 
obtaining of simple integral equations that are approximately 
quantitatively descriptive of boundary-layer phenomena 

Before a general treatment of the simplified equation, its appli- 
cation to the case of flow along a flat plate, where V is constant 
and the pressure term vanishes, will be briefly treated. The 
solution to the approximate equation for this case is of the form 


) . 


u 2 s , 
= (e~*) ds (33) 
i wee 


s = (y/2V2) V V/vx 


where 


In Fig. 1 this function is compared with the well-known exact 
solution for this case. For further comparison, the shear in- 


tensity at the surface, 
To) = pv(OUu/OY)y=0 


and the boundary layer-thickness, 


d* =f, {1 (u V)| dy 
0 


are compared with the exact solution in Table 1. 

In order to obtain the more general solution, the pressure term 
is written as a function of the velocity of the frictionless flow adja- 
cent to the boundary layer by means of Bernoulli’s equation 
gand Y = 


Further, the substitutions are made, u = 2V 4 


TABLE 1 


Shear Intensity Boundary-Layer 


(atx = a) Thickness, d* 
Exact solution 0.322p Vv V3/a 1.73 Vva/V 
Proposed approximation 0.38999 VvV3/a 1.6 Vva/V 
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yvV V/2p, giving 
0qg/Ox = 07¢/0¥? (4) 
It has already been shown that, for a constant velocity, V, 
u/V = f(y V V/va) 
Transforming this into the notation of Eq. (4), this becomes 
2V+q = Vf(¥vV2/a) 
Differentiating, 
(0g/OV)y=o = (0u/O¥V )y=9 = (0.899) VY 2/a 


For V is substituted (dV /dx)dx, the increment of velocity which 


arises at the point x. As the shear at the point x = 0 is desired, 


the relation a = (b — x) is substituted. Following this, sum 
mation of the incremental effects yields 
oq *b (dV/dx) 
: = ().399YV2 dx 
OV y= 0 V(b x) 


From this, it immediately follows that the shear at the point ) 
is given by the equation 
"6 (dV /d») 


dx (5) 


C, (0.899 pV Vv } 
0 Vv (b x) 


(to)x=b = 
where C; is a dimensionless coefficient to correct for the error of 


the original simplifying assumption. An expression for the 
boundary-layer thickness, obtained in a similar manner, is found 
to be 


. mg 
1.60 a , 


d* = (, Vv \ vd V(b x)(dV/dx) dx (6) 
where C; is a correction coefficient. The condition for separation 
is that, in Eq. (5), ro = 0. In both equations, the origin for x is 
the leading-edge stagnation point. 
constants so as to give exact results when the equations reduce to 
the case of a flat plate, C; = 0.806 and C, = 1.08. 

A second check on the validity of these equations is the calcu- 
lation of the separation point on a cylinder, assuming the pres 
sure distribution on its surface is that of ideal, symmetrical, fric 


If C; and Cy are evaluated as 


tionless flow. A series approach, based on the exact equation 
(Schlichting, p. 66+), yields a value of 92°. 


between 93° and 94°, measured from the leading stagnation point. 


Eq. (5) gives a value 


+ Schlichting, H., Lecture Series ‘‘ Boundary Layer Theory,’’ Part I—Lam 
N.A.C.A. T.M. No. 1217, April, 1949 


general reference, giving both sources and detailed discussion of the points of 


inar Flows, This paper serves as a 


boundary-layer theory referred to in this note 


Lateral Buckling of Beams in Bending and 
Compression 


F. DiMaggio, A. Gomza, W. E. Thomas, and M. G. Salvadori 
Department of Civil Engineering, Columbia University, New York 
April 21, 1952 


_— TOTAL ENERGY of an I beam buckled laterally under the 


action of normal loads and longitudinal thrusts is give. 


by! 3 
, 8 Ff Dr? [' H f' 
V= (u”)? dz 4+ (B")2 dz 4 (B’)? ds 4 
2 e o } e 0 2 e 0 
* P "S 
MBu" dz (u’)? ds (1) 
0 2 Jt 
where 
B = minimum flexural rigidity of beam 
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u = lateral deflection of beam in the plane of minimum rigid. 
ity 

/ = beam length 

D = flexural rigidity of one flange in its own plane 

h = beam depth 

8 = angular rotation of section of abscissa 

H = C—(PI,/A) = equivalent torsional rigidity of beam 

C = torsional rigidity of beam 

P = thrust, uniformly distributed on beam cross section 

I, = polar moment of inertia of beam cross section about its 
centroid 

A = cross-sectional area of beam 

M = bending moment in the plane of maximum flexyr, 


rigidity at section of abscissa 2 
= abscissa measured along the beam axis with origin 
one of its ends 
The cases considered in this note are among those for whict 
the bending moment diagram is linear, with values J/, at > = 
and M7, = rM.atz = 0, 
M(s) = Mz [r+ (1 r)s/l ) 
BEAMS UNDER 
Moments?* 


(A) Simmp_y Supported I UNEQUAL Enp 


The end sections of a ‘‘simply supported’’ I beam are free to 
rotate around their principal axes but are prevented from twisting 
around the longitudinal beam axis. Thus, the boundary condi 


tions for u and £ are, in this case, 


“(i) = w"(l) = 01 
Ail) = B*(l) = 0} 


“"(0) = 


B"(0) = 


u(O) = 
B(O) = 


and wv and 8 may be expanded into sine series, 


=> 


n 1 


a, sin naxx; B= > bn sin nxx (x = 2/l ! 
n 1 

In the absence of thrust, setting ? equal zero in Eq. (1), th 

equivalent torsional rigidity // becomes equal to the torsional 


The Ray 
taking into a 


rigidity C, and the last term of the equation vanishes 
leigh-Ritz method applied to Eq. (1) with ? = 0, 
count an increasing number of terms of the series (4), gives th 


results of Table 1. The critical values of the parameter, 


K = M/~V BC 
in terms of the parameter, 

I?/a? = 212?C/Dh? (6 
have been computed for r = 1, '/s, 0, 1/,, 1, taking into ac 


count 1, 2, 2, 3, and 4 terms of the series, respectively. The 


values of K corresponding to /?/a? = © (rectangular cross-sectiot 
beam) are within less than 1 per cent of the rigorous values com 
puted by means of Bessel functions® and listed in Table 1. Since, 
for /?/a? # 


the series taken into consideration is percentagewise very small 


o, the maximum change in A due to the last term ol 
it may be assumed that the values of Table 1 for r = 1, '/2, 0 
—1!/5, —1 have an accuracy of the order of 1 per cent. The re- 
maining values of K of Table 1 have been obtained by paraboli 
interpolation ;®° they have errors of the order of 2 per cent fot 
/?/a? = o, and it may be presumed that interpolation errors ot 
the same order of magnitude occur for the remaining values of 
[2 /q2,* 


(B) SrmpLy SuPPORTED, RECTANGULAR CROSS-SECTION BEAMS 
UNDER UNEQUAL END MOMENTS AND THRUSTS’ 


The total lateral buckling energy of this type of beam is ob 
tained by letting D = Oin Eq. (1). The boundary conditions for 


uand £ are, in this case, 


* Interaction curves for the lateral buckling of I beams under unequal en 


moments and thrusts will be published elsewhere in the near future 
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TABLE | 





- ; Ml _. ia 2/?¢ WV 
Values of K . Versus — = —— andr = 
B¢ a* h?D M, 
P/a 0.1 1.0 2.0 6.0 10 16 24 32 4() 100 
1.0 31.40 10.36 7.66 5.11 4.43 $4.00 3.43 3.59 3.51 3.29 3.14 
0.9 33.80 11.16 8.25 5.50 4.76 4.29 $.00 3.84 3.75 3.50 3.3 
0.8 35.75 11.81 8.73 5.82 5.03 4.54 $23 1.06 3.97 3.70 3.48 
0.7 37.53 12.40 9.16 6.11 5.28 ‘77 4.44 4.27 Bes 3.89 3.68 
0.6 39.35 13.01 9 61 6.41 5.55 5.01 1.67 41.49 4.38 $1.10 3.89 
0.5 $1.40 13.68 10.11 6.75 5.84 §.27 +. 92 4.73 $.62 1.33 t.12 
0.4 13.82 14.47 10.70 7.14 6.18 5.58 5.21 5.01 +.89 4.59 1.37 
0.3 16.69 15.41 11.39 7.60 6.58 5.93 5.54 5.33 5.20 4.87 +. 64 
0.2 50.07 16.52 12.21 8.14 7.04 6.34 5.92 5.69 5.55 5.20 4.93 
0.1 53.95 17.79 13.14 8.75 7.57 6.81 6.35 6.10 5.95 5.55 5.24 
0) 58.30 19.22 14.19 9.44 8.16 7.33 6.83 6.55 6.38 5.94 5.56 
ee 63.02 20.77 15.33 10.18 8.80 7.89 7.34 7.03 6.84 6.35 5.90 
-~(). 2 67.99 22.40 16.52 10.97 9.46 8.47 7.88 7.54 7.33 6.77 6.25 
—0.3 73.08 4 O6 17.74 11.76 10.14 9.07 8.41 8.04 7.82 7.20 6.61 
—().4 77.92 25.67 18.92 12.58 10.80 9.65 8.94 8.54 8.30 7.61 6.97 
—0.5 $2.39 7 . 14 20.00 13.24 11.40 10.18 9.42 9.00 8.74 8.00 7.32 
—().6 86.18 28.38 20.91 13.84 11.91 10.64 9.83 9.46 9.12 8.34 7.66 
-(0.7 88.79 29 .26 21.56 14.27 12.29 10.98 10.18 9.70 9.42 8.61 7.96 
~().8 89.97 29.66 21.86 14.48 12.48 11.16 10.30 9.88 9 59 8.78 8.19 
—().9 89.23 29.42 21.71 14.40 12.48 11.138 10.27 9.89 9.61 8.83 8.25 
—1.0 86.07 28.40 200. 98 13.95 12.07 10.85 10.01 9.69 9.43 8.72 8.03 
aaa a puted by Bessel functions® and listed in Table 2. Since the value 
Vol ‘ pp 7 of K for p = 1 equals zero for all values of r and the last term of 
Values of K = ; Versus p = and r = the series considered in the computations introduces a small 
BH 2 V 
V . re change in the value of K, it may be presumed that the values of 
Np () 0 2 0.4 0.6 08 1.0 K corresponding to r = 1, '/2, 0, : -] have errors of the 
10 3°14 2 8] 2 43 1.99 1.40 0 order of 1 per cent, while the interpolated values have errors of 
0.9 3.3 2.99 2.60 2.13 1.49 0 the order of 2 per cent, as shown by the case » = 0 when com- 
0.8 3.48 3.16 2.79 2.20 1.58 0 pared with the results of a rigorous solution.® 
0.4 3.68 3.33 2.89 2.3% 1.66 0 : ‘ = we Kv 
06 2 9g 3 5] 2 04 2 49 1.76 0 It is worth noticing that the curve K versus p for r = —1 pre 
3.8! 3.! ; 2.4$ f : sce as Bod aa : : : 
05 112 2 7) 2 9 2 63 1.86 0 sents a kink at p = 0.76, indicating a change in configuration for 
0.4 1.37 3.92 3.40 2.79 1.98 0 this value of the dimensionless thrust.* 
0.3 4.64 $.16 3.62 2.97 2.12 0 
0.2 1.95 + 44 3.88 3.19 2.28 0 
0.1 5.24 1.75 $.16 3.44 2.48 0 (C) CANTILEVER BEAMS OF RECTANGULAR CROSS-SECTION 
0) 5.56 5.08 1.48 3.73 2.70 0 UNDER THRUST AND SHEAR® 
—().] 5.90 5 1.83 4.05 2.96 0 
0.2 6.25 0 5.20 +. 40 3.25 0 In this case, indicating by s = 0 the free end of the beam and 
—0.3 6.6 ( 5.58 +.78 3.0 0 : 3 : . 
0 ‘ 4 ne a = oP 18 . 04 a by Q the shear applied at the centroid of the free section, Eq 
».Y4 ) oes I) ‘ Q o c 
-0.5 7.32 6 6.36 5.59 1.34 0 (2) gives 
—(0.6 7.66 7 6.73 6.00 4.78 0 ? 
-().7 7.96 7 7.06 6.41 5.26 0 M(s) = Qs = Moz/I (10) 
—().8 8.19 7 7.3 6.80 5.77 0 : . : : 
~1). 9 8.25 7 7 5A 7 14 6 31 0 Hence, Eq. (1) with D = 0 gives the total lateral buckling energy 
—-1.0 8.038 7 7.65 7.44 6.89 0 of the beam, provided r be set equal to zero in Eq, (2) 
The differential equations and boundary conditions satisfied 
F , by u and 8 may be derived from Eq. (1) by the calculus of vari 
u(0) = u"(O) = ul) = u"(1) = O7 = . 
(7) ations and are 
B(0) = B(l) = O \ 
gh Te Bu + Pu" M Bg)" =01 
ind wand 8 are expandable as in Eq. (4) : 7 . a ie (11) 
7 ‘ ; ae Mu 3, = 0 \ 
lable 2 contains the critical values of the parameter, ‘ 
K = Mil/V/BH eo = 
versus the parameter u(1) = u'(1) = u"(0) = BuO) + Pu(O) + QAO) = 0 | (12 
‘ , a 2) 
A(1) = B’(0) = 0 \ 
pb = P/P. = Pi?/x*B (9) 
; P . respectively. 
for? = 1, 1/30, —1/- 1, obtained by the Rayleigh-Ritz method P ‘sepals. 
9900 . The critical values of the parameter, 
ind computed by means of 1, 2, 2, 3, and 4 terms of the series, 
respectively. The other values of K in Table 2 were obtained by gq = O/Qer. = QI?/4.013 V BH (13) 
parabolic interpolation.6 The values of K for p = 0 have errors 
of less than 1 per cent when compared with the true values com * See footnote on p. 574 
TABLE 3 
Values of g = Q/?/4.013V BH Versus p = 4P/?/x*B 
p ] O.8 0.6 0.4 0.2 0 —(0).2 —().4 —1.0 
g 0 0.4776 0.6623 0.7984 0.9078 1.0000 1.0782 1.15038 1.3227 
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versus the parameter, 


p = P/P, = 4PP/xB (14) 
appearing in Table 3 were computed by elimination of 6 and ex- 
pansion of u into a power series. The curve q versus p is ap- 
proximated with errors of less than 2 per cent by the parabola, 
2.25 ~ 
p +q=1 (15) 
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Fluid Motions Whose Kinematics Are 
Independent of Compressibility of Fluid 


H. G. Loos 


National Aeronautical Research Institute, Amsterdam, Holland 
April 21, 1952 


ly A PAPER BY PARSONS,! it is investigated in what cases the 


streamline pattern of a compressible and an incompressible 
fluid will coincide. Parsons finds that this occurs if both the 
velocity and the vorticity are constant along stream lines. 

It is proved by Parsons that the only families of stream lines 
fulfilling these conditions are the concentric circles and the par- 
allel straight lines. 

This can be shown in a more simple way by transforming the 
Helmholtz equation in free stream-line coordinates. The cir- 
culation along the infinitesimal element ABCD (AB and CD 
along the parallel stream lines, BC and DA normal to it) is just 


(0/On) (qés)+ bn 
so the vorticity D 


(qés)* dn 
on 


bn - ds 
re) 


Og On 


On ra bs 


Og 60 Og 


+ q 
On 65s on R 


As q, 0g/On, and w are 
constant along the stream 
lines, it follows at once 
that the 
must be constant along the stream lines 


curvature 1/R 
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A Graphical Method for Determining Exact 
Boundary-Layer Temperatures 


J. J. Martin 
Aerophysics Laboratory, North American Aviation, Inc., 
April 11, 1952 
- IN DETERMINING boundary-layer temperatures, the spe- 

cific heat (c,) of air is considered constant, the boundary- 
layer temperature (7'g,) is given by 

Tat = To + 1,’ ( Vo?/2gJcp) (la) 

where the recovery factor n,’ = n,'(Tgr, T) and must be deter- 
mined for each (7%,, JT.) pair. Eq. (la) is familiarly written as 


Ta, = To {1 + 7,'[(y — 1)/2]M?} (1b) 


On the other hand, if variable specific heat (c,’) is allowed, then 
nr may be considered constant (for a given type flow) and 


*T RI eT « V, 9 
a dT = cp’ dT + 7, z (2a) 
0 2gJ 
*T aL V2 
or ed =_— (2b) 
T « 2¢J 


With the aid of the Gas Tables,* Fig. 1 was prepared with T in 
°R., and Vo in feet per second. In light of recent experiments 
showing 7,’ of the order of 0.6, values in the range 0.1 < , < 10 
were included for illustrative purpose. 

As an example, let 7. = 400°R. (approximately the tempef- 
ature in the isothermal region of the atmosphere), Vo = 4,000} 
ft. per sec. (= My) = 4.1 for T. = 400°R.), and 7, = 0.9 (cor 
responding, approximately, to turbulent flow); then 7s. = 
1,560°R. 1.40 and 9, = 1’, Eq 


(1b) gives Tpn = 1,552°R. 


On a similar basis, with y = 


* Keenan, J., and Kaye, J., Gas Tables, 2nd Ed., Table 1; John Wiley & 


Sons, Inc., New York, 1949 
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